NEW  YORK  UNIVEUSSTY, 
INSTlTf   :        .      ■    THEM ATICAL  SCIENCES 


'iS    Wa.eHV  n.LL,  Hlu    Y..I1   1|  M  Y' 


AEC  Computing  and 
Applied  Mathematics  Center 


AEC  RESEARCH  AND  DEVELOPMENT  REPORT 


PHYSICS  AND 
MATHEMATICS 


NYO-2538 


PAPERS   PRESENTED  AT   THE  SECOND 

INTERNATIOr^TAL  CONFERENCE  ON  THE 

PEACEFUL    USES  OF  ATOMIC  ENERGY. 

GENEVA,    SEPTEMBER  1958 

'.       January  30,.  1959 


Institute  of  Mathematical  Sciences 


NEW  YORK    UNIVERSITY 

NEW    YORK,    NEW    YORK 


{ 


!r* 


*  T^^^^i'i 


This  report  was  prepared  as  an  account  of  Government  sponsored  work.  Neither 
the  United  States,  nor  the  Commission,  nor  any  person  acting  on  behalf  of  the 
Commission: 


A.  Makes  any  warranty  or  representation,  express  or  implied,  with  respect  to 
the  accuracy,  completeness,  or  usefulness  of  the  information  contained  in 
this  report,  or  that  the  use  of  any  information,  apparatus,  method,  or 
process  disclosed  in  this  report  may  not  infringe  privately  owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to  the  use  of,  or  for  damages  resulting 
from  the  use  of  any  information,  apparatus,  method,  or  process  disclosed 
in  this  report. 

As  used  in  the  above,  "person  acting  on  behalf  of  the  Commission"  includes 
any  employee  or  contractor  of  the  Commission  to  the  extent  that  such  employee 
or  contractor  prepares,  handles  or  distributes,  or  provides  access  to,  any  in- 
formation pursuant  to  his  employment  or  contract  with  the  Commission. 


UNCLASSIFIED 


AEC  Computing  and  Applied  Mathematics  Center 

Institute  of  Mathematical  Sciences 

New  York  University 


PHYSICS  AND  NYO-2538 

MATHEMATICS 


PAPERS   PRESENTED  AT   THE  SECOND 

INTERNATIONAL  CONFERENCE  ON  THE 

PEACEFUL    USES   OF  ATOMIC  ENERGY, 

GENEVA.    SEPTEMBER  1958 


January  30,    1959 


Contract  No.    AT(30-1)-1480 


1  - 


UNCLASSIFIED 


NYO-2538 

PHYSICS  AND 

MATHEMATICS 


PREFACE 

This  report  contains  corrected  manuscripts  of  five  papers  sub- 
mitted by  the  New  York  University  Sherwood  group  to  the  Second  Inter- 
national Conference  on  the  Peaceful  Uses  of  Atoinic  Energy  in  Geneva, 
September  1958.     Also,    some  references  to  previously  classified  litera- 
ture have  been  added. 


Harold  Grad 

January  14,    1959 


/ 


NYO-2538 

TABLE  OF  CONTENTS 

Page 
Preface 2 

Hydromagnetic  Shock  Waves  in  High- Temperature  Plasmas  .....  4 

C.    S     Gardner,    H.    Goertzel,    H.    Grad,    C.    S,    Morawetz, 
M.    H.    Rose,    H.    Rubin 

Problems  in  Magnetohydrodynamic  Stability 29 

J.   Berkowitz,    H.    Grad,   H.    Rubin 

Hydromagnetic  Equilibria  and  Force- Free  Fields 68 

H.    Grad  and  H.    Rubin 

Thermonuclear  Reaction  Rate  in  a  Discharge 94 

H.   Grad 

Cusped  Geometries 108 

J.    Berkowitz,    K,.    O,    Friedrichs^    H,    Goertzel,    H..    Grad, 
J.   Killeen,    E.    Rubin 


3  - 


NYO-2538 


HYDROMAGNETIC  SHOCK  WAVES  IN  HIGH- TEMPERATURE  PLASMAS 

C.    S.    Gardner,    H.    Goertzel,    H.    Grad, 
C.    S.    Morawetz,    M.    H.    Rose,    H.    Rubin 

1.     INTRODUCTION 

Shock  waves  are  one  of  the  most  important  tools  for  heating  ordinary 
gases  and  for  creating  a  plasma.     It  is  still  an  open  question  whether  shock 
heating  will  have  as  vital  a  position  for  high  temperature  thermonuclear  plas- 
mas.    The  theoretical  results  are  not  quite  complete,    and  experimental  con- 
firmation may  take  a  while.     However,    the  present  results,    summarized  m 
this  paper,    do  justify  a  considerable  degree  of  optimism  as  to  the  suitability 
of  a  certain  range  of  shocks  for  heating  to  thermonuclear  temperatures. 

In  classical  nondissipative  gas  dynamics,    it  is  observed  that  a  finite- 
amplitude  compressive  wavefront  progressively  steepens  until  it  becomes  mul- 
tivalued and  thereby  physically  meaningless.     This  mathematical  catastrophe 
can  be  avoided  by  the  introduction  of  dissipation,    usually  as  heat  flow  and  vis- 
cosity.    The  widening  effect  of  the  dissipative  mechanism  becomes  more  pro- 
nounced as  the  wavefront  steepens,    and  a  wave  shape  of  permanent  form  is 
reached  as  a  balance  between  the  effects  of  nonlinearity  and  dissipation.     The 
re^^rence  length  for  this  transition  zone  is  the  mean-free- path,    and  its  actual 


thickness  ranges  from  one  to  several  mean-free-paths  except  in  the  case  of 
very  weak  shocks  which  are  wider.  This  shock  layer  is  therefore  experi- 
mentally indistinguishable  from  a  discontinuous  wave  front  in  all  conventional 


Grad,    H.  ,    The  Profile  of  a  Steady  Plane  Shock  Wave,    Communications  on 
Pure  and  Applied  Mathematics,    5:  257-300,    (1952)  and  additional  references 
quoted  there. 
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gas  dynamical  applications.     The  state  on  one  side  of  the  shock  can  be  com- 
puted from  the  state  on  the  other  side  (and  one  additional  parameter,,    the 
"strength"  of  the  shock)  by  direct  use  of  the  laws  of  conservation  of  mass, 
momentum,    and  energy.      The  particular  nature  of  the  dissipative  mechanism 
determines  the  shape  of  the  transition  zone  but  not  the  end  states.     This  for- 
tunate circumstance  makes  it  possible  to  solve  gas  dynamical  flow  problems 
with  the  relatively  simple  nondissipative  fluid  equations  by  piecing  together 
smooth  solutions  at  appropriately  selected  shock  fronts. 

In  a  conducting  fluid  (plasma),    the  results     obtained  by  conventional 
theory  are  essentially  the  same.     In  particular,    inclusion  of  a  third  dissipa- 
tive mechanism,    finite  electrical  conductivity     leaves  the  thickness  on  the 
order  of  the  mean-free-path  or  larger.     The  presence  of  a  magnetic  field 
complicates  matters  considerably,    but  the  essential  features  are  unaltered. 
In  high-temperature  plasmas,    the  mean-free-path  can  be  many  times  larger 
than  the  size  of  the  apparatus.     The  apparent  conclusion  is  that  one  would 
not  be  able  to  observe  a  shock  wave  in  a  thermonuclear  plasma  of  laboratory 
scale.     This  would  be  unfortunate  on  two  counts       Mathematically,    it  would 
necessitate  the  use  of  complicated  differential  equations  to  solve  problems 
that  could  otherwise  be  solved  explicitly  by  algebraic  means.      Physically,    it 
would  seem  to  eliminate  one  of  the  most  promising  means  of  irreversible 

heating  of  a  plasma.     Fortunately,    the  conventional  theory,    i   e.    the  conven- 

2 
tional  modification  of  the  stresses  and  heat  flow  in  a  magnetic  field"    is 

completely  inapplicable  when  the  gyro  radius  is  stnaller  than  the  mean- free- 
path. 


Marshall,    W.  ,    The  Structure  of  Magneto- Hydrodynamic  Shock  Waves,    Pro- 
ceedings of  the  RoyaFSocietj^lLon^^  (1955). 

2 

Chapman,    S.    and  Cowling,    T,    G.  ,    The  Mathematical  Theory  of  Non-uniform 

Gases,    Cambridge  University  Press,    p     338,    (1939). 
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With  the  mean-free-path  eliminated,    we  naturally  turn  to  the  Debye 
length  or  the  gyro  radius  as  a  suitable  reference  length.         The  significant 
length  parameter  turns  out  to  be  neither  of  these;     it  is  the  speed  of  light 
divided  by  the  plasma  frequency.     This  is  roughly  intermediate  between  the 
electron  and  ion  gyro  radii.     The  problem  is  to  exhibit  an  irreversible  mech- 
anism which  operates  on  this  scale  of  length, 

To  isolate  the  essentials  of  the  problem,    we  take  the  mean-free- 
path  to  be  infinite.     An  electron  or  ion  feels  only  the  electromagnetic  field 
of  the  averaged  distant  charges  which  are  represented  by  the  source  terms 
in  Maxwell's  equations.     This  is  a  so-called  self- consistent  field  formula- 
tion.    The  crucial  point  to  realize  is  that  irreversibility  has  not  been  elim- 
inated by  removing  the  collision  term  in  the  Boltzmann  equation,    it  has 
merely  been  concealed.     We  can  still  rely  on  the  fundamental  irreversible 
mechanism  described  heuristically  by  Gibbs.     A  special  case  in  the  applica- 
tion of  this  mechanism  to  plasmas  has  been  referred  to  as  Landau  damping. 
The  problem  we  consider  is  more  subtle  and  the  solution  has  been  frequently 
held  to  be  undamped.     The  basic  property  which  distinguishes  every  irrevers- 
ible  mechanism  is  that  two  microscopic  initial  states  which  are  so  close  as 
to  be  indistinguishable  may  eventually  emerge  very  far  apart.     Short-range 
intermolecular  collisions  (which  we  have  chosen  to  igiiore)  exhibit  this  prop- 
erty very  clearly.     This  mechanism,    sometimes  called  phase  mixing,    is 
intuitively  simple  but  frequently  offers  great  mathematical  difficulties. 


To  be  sure,    the  Debye  length  already  exists  as  an  alternative  size  even  with- 
out the  magnetic  field.     However,    it  can  be  shown  that  a  compressive  motion, 
e.  g.    one  produced  by  a  piston,    will  widen  indefinitely  m  the  absence  of  a  mag- 
netic field  until  particle  collisions  take  over.      Paradoxically,    collisions,    which 
are  the  source  of  the  conventional  dissipative  widening  effect,    are  required  to 
effect  a  state  in  which  gas  dynamics  becomes  applicable,    i.  e.  ,    in  which  com- 
pression waves  steepen. 


The  practical  question  which  we  must  answer  is  whether  a  violent 
irreversible  change  of  state  can  be  produced  within  a  thin  layer  of  some 
permanence.     The  theoretical  question  is:  does  the  shock  exist?     In  addi- 
tion,  the  heating  of  ions  and  electrons  separately  must  be  found.    It  is  nec- 
essary to  solve  the  shock  transition  problem  in  a  high-temperature  plasma 
not  only  to  establish  its  observability  (i.  e.    its  thickness)  but  even  to  find 
the  correct  jump  conditions.     In  the  absence  of  collisions,    there  is  no  rea- 
son for  the  state  after  the  shock  to  be  in  thermodynamic  equilibrium,     spe- 
cifically,   the  ion  and  electron  temperatures  can  be  different.     Conservation 
of  mass,    momentum,    and  energy  can  only  predict  the  mean  temperature. 
We  find  that  it  is  possible  to  choose  conditions  such  that  the  shock  will  heat 
the  ion  gas  much  more  than  the  electron  gas. 

In  the  following  sections,    a  number  of  approximate  two-fluid 
theories  are  described,    also  several  approximations  to  a  complete  self- 
consistent  field  solution,    as  well  as  combinations  of  fluid  and  particle 
analyses.     The  cumulative  evidence  seems  to  be  quite  strong  for  the  exist- 
ence of  shock  waves  with  the  properties  described. 

2.     FORMULATION  OF  THE  PROBLEM 

We  consider  a  steady  flow  in  the  positive  x-direction  in  which  a 
shock  transition  occurs  from  a  constant  state  at   x  =  -  oo    (front  of  the  shock) 
to  another  constant  state  at   x  =  +  oo    (back  of  the  shock).     All  macroscopic 
quantities  are  assumed  to  be  functions  of  x   alone.     A  magnetic  field,     B, 
points  in  the  z- direction.     It  is  accompanied  by  a  current,    J,    in  the  y- 
direction 


We  use  rationalized  MKS  units. 


-  7 


dx 


(1) 


The  electric  field   E   has  a  component    E      in  the  y-direction  and  possibly  a 


yf 


V.Ey.j 


u.E, 


■♦    X 


B 


Fi9.l 


charge  neutrality  for  all  x: 


component    E      (charge  separation)  in  the  x- 
direction.     From   aB/9t  «  -  curl  E  »  0,    we  con- 
clude that    E      is  a  constant,    independent  of  x. 

y 

Velocity  components  are   (u,  v)    in  the   (x,  y)    direc- 
tions;   in  a  two-fluid  theory,   the  current,    J, 
arises  from  the  velocity  components    v.      We  dis- 
tinguish between  ion  and  electron  quantities  by 
subscripts,    u   ,    u   ,    etc.     At    x  «  +  oc  and   x  »  -  co 

we  have   v     *  v     «  J  =  0   and   E     «  0.      We  assume 

+  -  X 


n       »    n       «    n 


(2) 


(n  is  the  number  density).     This  is  valid  if  the  Debye  length  is  small  compared 
to  the  scale  of  the  shock  transition.     For  most  of  the  arguments  in  this  paper, 
this  is  a  convenience  rather  than  a  necessity.     From  the  two  laws  of  conserva- 
tion of  mass  (or  number),    we  have 


n   u       *    constant 


n  u       *    constant      ,    (3) 


whence,    for  all  x,    we  conclude  that 


u       «    u       =    u 

+ 


(4) 


At  both   X  «  +  00   and  x  *  -  oo    we  must  have 


E       »    uB 

y 


(5) 
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The  state  in  front  is  assumed  to  be  not  only  constant  but  in  thermo- 
dynamic equilibrium,     T     ■-  T    .      At   x  =  +  oo  we  may  have  unequal  tempera- 
tures,    T     ji  T    .      Of  course,    collisions  would  eventually  thermalize  the  two 
gases,,    but  this  is  not  part  of  the  problem  considered. 

We  introduce  the  following  self-explanatory  notation; 

p.     =    nkT.     =    p.R.T.      ,  p.     =    m.n      ,       R.     =    k/m.     ;         (6) 

1  1111  11  1  1 

P     -    P^   +    P_     ,         p    ^   p^   +   p_     ,         T    =  |(T^   +   T_)     , 


p    =    pRT      ,  R    ^   k/m    ,        m    =  "^("^+  +    m   )        ; 

pu=pu+pu        ,        pv   '^   p    V     +  p    V      ,        J   =  ne(v,    -    v   ) 


(7) 


(8) 


3.     JUMP  CONDITIONS.     HEATING 

From  a  particle  viewpoint,   the  z- component  of  velocity  of  each 
particle  is  a  constant  of  the  motion;     only  the  x-  and  y- components  matter. 
Hence  any  fluid  theory  which  describes  the  situation  must  be  strictly  two- 
dimensional,    i.  e.   the  internal  energy  per  unit  mass  is  given  by   e  =  RT 
rather  than    3/2  RT,    and  the  gas  constant  (ratio  of  specific  heats)  is  7  =  2. 
The  laws  of  conservation  of  momentum  and  energy  —  mass  has  already  been 
taken  care  of  by  equation  (3)  —  are  expressible  as 

[pu"    +    p   +   B"/2m]     =    0  .  [pu(e   +  ^u    »    +   pu   +    E   B/fu]     =    0     ;  (9) 

the  sytubol   [Q]  =  Q    -  Q      is  the  jump  m  Q  across  the  shock.     In  terms  of  the 
"total  pressure"  and  "total  energy" 

p*     =    p  +  B^/2m  ,  e*     =    e  +  B^/2iJLp      ,  (10;» 


the  jump  conditions  (9)  are  identical  to  those  with  no  magnetic  field, 

[pu"   +    p*]     =    0  ,  [pu{e*    +  |u^^    +   p*u]     =    0  (11) 

We  also  have  the  corresponding  equations  of  state, 

p    =    pe  ,  P'     ==    pe''"      .  (12) 

Thus,    in  terms  of  the  reduced  variables   p,    u,    p''',    e*  (and   T'''  =  p''' / P  if  desired; 
the  jump  conditions  are  the  saine  as  those  for  an  ordinary  gas  with  7=2.         To 
find  the  jump  in  B  and  in  the  original  variables    p,    e,    and  T   it  is  only  necessary 
to  supplement  the  reduced  jump  conditions  with  the  relation   [Bu]  =  0    (cf,    equa- 
tion (5)1. 

We  introduce  three  speeds'     the  sound  speed  a,    the  Alfv€?'n  speed  A, 
and  the  reduced  sound  speed  a"",     with  three  corresponding  Mach  numbers, 

a"     =    2RT      ,  M       -^    u/a 

a 

A^     =    B^/mp      ,  M^    -    u/A  (13,» 

(a*:-^     =    a^   +  A'^      ,  M*     =    u/a*      . 

We  recall      that  the  value  of  any  one  of  the  shock  strength  parameters    M''', 

Pjp    ,    u    /u, ,    p'//p*.    T''7T'''    fixes  all  the  others,      (Henceforth  we  shall  write 
1      o       o      1       1      o'       1       o 

M''~  instead  of  M''\  )  In  other  words,  there  is  only  one  dimensionless  param- 
eter characteristic  of  a  gas  shock.  On  the  other  hand,  there  are  two  dimen- 
sionless parameters  for  the  magnetic  shock.      These  can  be  chosen  as  M'''  and 

P.. 
3        =    -ir^^       ,  (14) 

°  b2/2m 

o 


For  example,    see  Courant,    R.    and  Friedrichs,    K.    O.  ,    Supersonic  Flow  and 
Shock  Waves,    Interscience  Publishers,    New  York,    (1948). 
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the  ratio  of  gas  to  magnetic  pressure  in  front.     Or,    we  can  take  two  of  the 
Mach  numbers  (13)  to  describe  the  shock. 


Although  the  pressure  and  temperature  ratios  across  a  shock,     p''7p''' 

and   T'/T   ,    can  become  arbitrarily  large  for  strong  shocks,   the  density  and 

velocity  ratios,    p  Ip      and  u    /u       are  at  most  three  for  a  gas  with  7=2. 

The  same  bound  therefore  holds  for  the  magnetic  field  ratio   B  /B     =  u    /u, . 

^  1       o         o      1 

2      2  2       2 

From  the  easily  verified  fact  that    pf/p*  >  u    /u      =  B    /B      and  from  the  rela- 

1      o         o      1  1       o 


111    mc    ccLoii^     vciiiic^    icn^^L    Liicit     p'/>^ 

tion 


Pi    ^1'  ^-^'^l 

—    =   — ^ (15) 

p  B  o 

o  o 

we  see  that   p, /p    >  pf/p*.      Moreover,    if  we  keep  the  strength   M      (or  B, /B   ) 
1     o        1     o  1      o 

fixed  and  let  /3      approach  zero,    we  obtain  a  limiting  value  for  |3     which  is  dif- 
ferent from  zero.     In  other  words,    if  the  magnetic  field  dominates  in  front  of 
the  shock,   then  even  a  shock  of  moderate  strength  can  heat  the  gas  enormously: 
T  /T     -^  00  ,      The  gas  is  heated  to  reach  an  energy  per  unit  volume  which  is  a 
certain  fraction  of  the  magnetic  energy,    independent  of  the  original  gas  energy. 
For  example,    a  shock  of  medium  strength,    M'''  =  2,    impinging  on  a  plasma  at 

"zero"  temperature,    jS     =0,    will  heat  it  to  a  value  (3,  =  1/4. 

o  1 

In  order  to  separate  the  ion  and  electron  energies  at  the  back  of  the 

shock  it  is  necessary  to  solve  the  flow  problem  for  the  transition  or  else  find 

an  additional  jump  condition.     We  shall  find  that  under  certain  circumstances 

this  extra  condition  is  that  the  electron  gas  is  approximately  adiabatic.     The 

word  adiabatic  is  used  in  two  different  but  equivalent  senses.     If  the  individual 

electron  motion  is  adiabatic  in  the  sense  that  its  magnetic  moment  is  approxi- 

1  0 

mateiy  constant,    we  have  —  mv   /B  =  constant.     For  each  electron  we  have 

^  2 

energy~B~p;    therefore,    for  the  whole  electron  gas  we  have   RT~p   or   p~p    , 


11 


which  IS  the  adiabatic  gas  law  for  7  -  2.      We  shall  discuss  the  justification 
of  this  adiabatic  assumption  later.     The  result  is  that  the  electron  tempera- 
ture  T      can  increase  by  at  most  a  factor  three;    the  ions  are  always  heated 

more,    and  if  (3      is  small,    the  ion  temperature   T      can  increase  very  much 
o  + 

more.     Specifically,    writing  t    for  the  ratio  of  the  electron  temperature  be- 
hind the  shock  to  that  ahead  (so  that    1  <  '^  <  3),    we  find 

+  0  O 


4.     ADIABATIC  TWO- FLUID  THEORY 

This  is  perhaps  the  simplest  theory  that  can  be  proposed  which  dis- 
tinguishes between  the  ion  gas  and  the  electron  gas.     They  are  each  individu- 
ally assumed  to  obey  conventional  nondissipative  fluid  equations  supplemented 
by  a  "Lorentz"  force.     We  know  beforehand  that  we  shall  not  find  a  solution 
which  joins  two  end  states  compatible  with  the  shock  conditions  since  the  latter 
imply  an  entropy  increase.     Nevertheless,    the  results  are  very  illuminating. 

The  meaning  of  the  term  "nondissipative"  gas  is  somewhat  flexible. 
One  can  assume  that  there  is  a  scalar  pressure  which  satisfies  the  adiabatic 

law 

9 
p    =    ap"      .  (17) 

Or  one  can  take  an  energy  conservation  equation  of  conventional  form  involv- 
ing a  scalar  pressure  and  no  heat  flow.     In  each  case,    there  is  some  ambiguity 
m  the  term  scalar  pressure  depending  on  whether  the  frame  of  reference  is 
the  speed  of  the  individual  gas  or  of  the  gas  mixture;     for  simplicity  we  adopt 
the  former  definition.     In  the  case  of  a  simple  gas,    the  energy  equation  has  the 
adiabatic  relation  (17)  as  an  integral.     For  the  mixture  of  two  gases,    the  t) 


two 
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postulates  are  distinct.     The  equation  of  energy  conservation  would  seem  to 
be  more  fundamental.     The  results  obtained  from  the  two  theories  are  quali- 
tatively similar.     Since  only  the  qualitative  features  of  the  solution  have  any 
bearing  on  the  shock  problem  as  posed,    we  shall  describe  here  the  analytically 
simple  adiabatic  theory. 

The  total  (i.  e.    sum  of  ion  and  electron)  mass  and  x-momentum  equa- 
tions can  be  integrated. 

pu  =  CT      ,  pu^   +   p  +  B^/2m     =    P      .  (18) 

Together  with  the  adiabatic  law  (17),    we  find  the  relation  between  B  and  u: 

2 
B^     «    2/u[P  -  au  -  ^]      .  (19) 

u 

The  difference  between  the  two  y-nriomentum  equations. 

4-  (p.uv  )     ■■-    ne.(E     -  uB)  (20) 

dx     1      1  1     y 

can  be  combined  with  the  Maxwell  equation  (1)  to  give 

d^B                2,1           1    ,,„  ^    ,    , 

— -  -    june    ( + '(B  -  E    /u; 

,2                      mm  y 

dx                            +           .  ^ 


2  dj 

^(B  -  E    /u)    =    --- 
u  y  dB 


(21) 


We  recognize  the  dimensional  factor  m  the  first  line  to  be  the  squared  ratio 

2 
of  plasma  frequency  to  the  speed  of  light;    7^      is  a  constant.     Since  u  is  known 

implicitly  through  equation  (19)  as  a  function  of  B,    equation  (21)  can  be  inte- 
grated.    The  behavior  of  the  solutions  are  immediately  evident  when  we  inter- 
pret  6(B)    as  a  potential  for  the  "motion"    B(x' .     From  Fig.    2,    we  see  that 
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there  is  a  family  of  periodic  motions  of  increasing  amplitude  around  the  poten- 
tial minimum  (back  of  the  "shock")  culminating  in  a  pulse  or  solitary  wave 


B 


Fig.  2 


Fig. 3 


which  starts  at  the  front  and  eventually  returns.     These  solutions  are  shown  in 
the  phase  plane.    Fig.    3,    and  are  sketched  as  functions  of  x   in  Fig.    4. 


Fig.  4 

The  entropy  can  be  chosen  as  the  parameter  which  distinguishes 
these  curves.     It  should  certainly  be  expected  that  a  dissipative  mechanism 
would  produce  the  damped  oscillation  indicated  in  Fig.    4  by  the  dashed  lines, 
producing  a  legitimate  shock  transition  from  front  state  to  back  state.     This 
has  actually  been  verified  by  inserting  a  friction  term  in  the  momentum  equa- 
tions (ohmic  dissipation)  proportional  to  the  relative  velocity,   but  the  damping 
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length  is,    of  course,    the  mean-free-path. 

The  reference  length  for  the  periodic  oscillations  is  essentially 

h       ^ .  (22) 

2 
lune 

The  phase  length  or  distance  travelled  by  a  particle  in  one  gyro  period  is 

L      ^    m   u    /eB      .  (23) 

1  1     o 

2 
We  find  that   L    /L  -  eM,    and  L    /L  ^  M„  /e    where   e     =  m    /m    .      For  M  .    not 

A  +  A  -        +  A 

too  large  (also  not  too  small  -  for  weak  shocks,    the  actual  wavelength  becomes 
large  compared  to  the  reference  length,    L),    we  conclude  that  the  individual 
electron  orbits  are  quite  accurately  adiabatic    (L     <<  L),    whereas  the  ion  or- 
bits definitely  are  not.     Provided  that  there  does  exist  a  shock  transition  ap- 
proximated by  this  theory,    the  electron  gas  will  be  adiabatic;    the  ion  gas  will 
be  heated  more  and  in  some  cases  much  more.     For  the  thermonuclear  appli- 
cation,   the  most  interesting  shocks  would  seem  to  be  those  of  medium  strength. 

The  value  of  this  theory  is  that  it  strongly  suggests  that  we  look  for 
an  oscillatory  shock  transition  rather  than  the  classical  monotone  transition. 

More  elaborate  fluid  theories,    including  non-scalar  stresses  and  heat 
flow  have  been  attempted,    but  the  results  are  not  yet  definitive. 

5.     SELF- CONSISTENT  FORMULATION.     THE  DAMPING  MECHANISM 

If  the  fields    E(x)    and  B(x)    are  known,    the  motion  of  an  individual  par- 
ticle IS  governed  by  the  equations 


du  p  dx 

—     -    -(E   (x)    +   vB(x))      ,  u    -  -5^ 

dt  m        x  dt 


^     ^    ^(E     -uB(x))      ,  v  =  ^     . 

dt  m       y  dt 
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(24) 


(There  should  be  no  confusion  in  using  u  and   v  for  particle  velocities  as  well 
as  for  gas  velocities    )    The  second  equation  can  be  integrated  once  in  terms 
of  the  "vector"  potential. 


A(x)     ^ 

(25) 

v    ^-    —  [E    t  -  A(x)]  +  constant, 
m      y 

If  the  energy  of  the  particle  is  bounded     we  conclude  that  its  position  is  never 
more  than  a  bounded  distance  away  from  the  solution  of  the  implicit  relation 
E  t  -  A(x)  -  0      which  is  equivalent  to   dx/dt  -  E    /B.      In  particular,    all  parti- 

y  y 

cles  eventually  move  on  to  x  =  +  co 

In  principle,    equations  (24)  can  be  solved  for  any  given  coefficients 
E(x)    and  B(x'    and  for  any  given  initial  values       Now  suppose  that  there  is  a 
Maxwellian  velocity  distribution  of  particles  for   x   very  far  to  the  left,     A  solu- 
tion of  Liouville's  equation,    which  is  equivalent  to  a  solution  of  the  particle 
trajectories  for  all  initial  states,    yields  the  particle  distribution  for  all   x      By 
solving  this  equation  for  both  the  electron  gas  and  the  ion  gas  (with  the  given 
E(x)    and  B(x)),    we  can  compute  the  densities,    velocities,    and  current  as  func- 
tions of  X.      If  the  current  computed  in  this  way  is  consistent  with  the  original 
B(x)      and  if  the  ion  and  electron  densities  turn  out  to  be  equal  for  all  x,    then 
we  have  obtained  a  self- consistent  solution  to  the  problem. 

Any  verification  of  the  existence  of  a  solution  to  the  shock  problem 
just    posed    must    be    quite    subtle  There  is  a  very  simple  and  plausible 

argument,    based  on  Liouville's  theorem,    which  apparently  proves  the  non- 
existence of  a  shock  under  the  circumstances  considered,    except  that  it  de- 
pends on  certain  unwarranted  assumptions       Consider  the  molecular 
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Fig.  5 


distribution  function,    f(x,  y,  u,  v).      It  satisfies  Liouville's  equation  which 
states  that   f  is  a  constant  on  every  particle  trajectory.     At   x  »  -  co,    f  is 
Maxwellian;    in  particular,   the  lines  of  constant    f  are  circles  in  the  u,  v 
plane  (Fig.    5).     The  existence  of  a  constant  state  at   x  »  +  oo    implies  that   f 

is  isotropic  in  velocity;    the  lines  of  con- 
stant  f  are  again  circles,    but  centered  about 

u     rather  than  u    .      Now,    for  some  value  of 
1  o 

X  far  to  the  left,    consider  an  element  of  the 

four- dimensional  phase  space:    specifically, 

consider  the  product  of  an  area  element    dA 

in  physical  space  and  a  region   27rrdr  between 

two  circles  in  velocity  space;    f  is  essentially 

constant  on  this  domain.     Although  the  area  element   dA   may  be  deformed  in  a 

very  complicated  way  in  moving  to  x  «  +  oo,   the  projection  onto  velocity  space 

must  again  be  circular  at   x  »  +  oo    since  it  must  lie  on  a  curve  of  constant   f. 

The  only  question  is  which  circle  maps  onto  which.     According  to  Liouville's 

theorem,   the  four- dimensional  volume  remains  constant  during  the  motion. 

Just  as  the  overall  conservation  relation  n  u     «  n,  u     is  obtained,   a  similar 

o   o         11 

one  can  be  obtained  for  any  specially  identified  group  of  particles.     We  con- 
clude that  the  area  element    dA   is  compressed  by  a  factor  n  /n      in  moving 

1      o 

from   x*-oo    tox««+oo.      The  area  element  in  velocity  space  is  expanded  by 
the  same  factor.     This  implies  that  the  radius  of  the  original  circle  is  ex- 
panded by  a  factor    i  njn     ,    from  which  we  conclude  that  the  distribution  at 

1      o 

X  *  +  00    is  also  Maxwellian  at  a  temperature  which  is  higher  by  the  factor 

"i /n    .      But  this  is  an  adiabatic  compression  and  is  inconsistent  with  the  shock 

jump  conditions  unless  the  two  end  states  are  identical. 

This  argument  is  fallacious  in  exactly  the  same  sense  as  the  conven- 
tional arguments  which  pretend  to  demonstrate  that  no  irreversibility  can 
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appear  in  a  Hamiltonian  dynamical  system.     It  has  been  tacitly  assumed  that 
there  is  pointwise  convergence  of  the  distribution  function   f  to  its  assumed 
limit  at   X  »  +  00.      A  little  consideration  will  show  that  not  only  is  the  area 
element   dA   deformed,    but  the  image  of  the  circle  in  velocity  space  becomes 
extremely  deformed  (cf.    Fig.    6)  and  only  converges  weakly  to  its  limit;    two 

points  in  phase  space  which  are  very  close  at 
X  «  -  CO    can  become  widely  divergent  as   x  —  +  oo. 
A  circle  at   x  »  -  oo    becomes  smeared  onto  a  two- 
dimensional  region  at   x  «  +  oo. 


Fig.  6 


It  is  possible  to  generalize  the  concept  of 
constant  state  at  the  back  of  the  shock  in  the  follow- 
ing way.     If  the  ion  and  electron  distributions 
are  symmetric  about  the  mean  velocity  (even  though  not  isotropic),   there  will 
be  no  net  current.     Consequently,   there  exist  exact  self-consistent  solutions 
in  which  a  non-constant  microscopic  state  is  coupled  with  a  constant  macro- 
scopic state.     One  might  think  that  such  non-isotropic  states  are  equally  possi- 
ble as  limiting  states  behind  the  shock.     However,  it  is  plausible,  from  the  same 
qualitative  arguments  as  above,    that  such  states  are  unstable  unless  they  are 
isotropic. 

We  have  already  observed  that  the  case  3     «  0    (zero  temperature  in 

o 

front  of  the  shock)  is  particularly  interesting  from  the  point  of  view  of  heating. 
At  zero  temperature,    there  is  no  dispersion  of  ion  velocities,    and  only  a  single 
orbit  has  to  be  calculated.     If,    in  "the  particle  equation  (24)  we  set    d/dt  «  u(d/dx), 
we  observe  that  the  equation  of  motion  of  a  single  particle  is  identical  with  the 
fluid  equation  obtained  by  setting  p  »  0.      The  exact  self- consistent  field  prob- 
lem has  therefore  been  solved  in  this  special  case  in  Section  4.         There  exists 


This  zero-temperature  solution  has  also  been  obtained  by  Davis,    Lust,    and 
Schluter,   also  by  Mittelman,    and  probably  by  many  others. 
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no  self- consistent  shock  solution  (at  least  for  shocks  of  moderate  strengt,h; 
for  strong  shocks,  the  question  is  still  open        It  is  clear     a  priori,    that  we 
could  expect  no  shock  m  this  special  case      The  damping  tnechanism  which 
we  have  just  discussed  requires  a  dispersion  of  velocities  about  the  mean, 
1    e.    a  finite  temperature.     We  do  derive  the  useful  piece  of  information  that 
the  shock  width  (i    e     the  damping  distance)  will  be  much  larger  than  the  ref- 
erence length,     L.     (equation  (22);  if  /3     is  very  small. 

In  the  remainder  of  this  paper  we  shall  concentrate  mainly  on  medi- 
um strength  shocks  for  which  the  electron  motion  is  assumed  to  be  adiabatic 
This  converts  the  problem  to  a  combination  of  one  essentially  fluid  equation 
(for  the  electrons^  and  one  particle  equation  (for  the  ions».     One  particularly 
simple  form  of  the  adiabatic  electron  relations       is     for  exainple 

u(xt     --    u   (x)     =-    E    /B(x)      ,  V   (x.     ^    -  E   (xWB(x) 

y  ■  -  X 

(26) 

n(x>     -    constant /u(x,f      ,  J       =    -  nev 

The  x- component  of  velocity  and  the  electron  contribution  to  the  current 
(which,    under  certain  circumstances,    can  be  taken  to  be  the  entire  current) 
are  explicitly  given  in  terms  of   E    (x»    and   B(x)        To  this  approximation,    the 
magnetic  field  self- consistency  relation  becomes 

r 

B(xl     =    constant    \   E   (x)dx 

X 


This  is  the  zero-order  guiding  center  approximation  to  the  current,     see 
Alfven,    H   ,    Cosmical  Electrodynamics,    Clarendon  Press,    Oxford     (1950) 
Kruskal     M     D    ,    The  Gyration  of  a  Charged  Particle,    NYO-7903,    (1958«  and 
Berkowitz     J    ,    and  Gardner,    C     S       On  the  Asymptotic  Series  Expansion  of 
the  Motion  of  a  Charged  Particle  in  Slowly  Varying  Fields,    NYO-7975,    (1957), 
Equation  (26)  can  easily  be  improved. 
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The  problem  is  to  find  an   E   (x)    such  that  the  resulting  ion  density  agrees 
'^  X 

with  the  above  computed  electron  density. 


6.     EXPANSION  IN  THE  MASS  RATIO       REBOUNDS  AS  A 
MECHANISM  FOR  DAMPING 

The  results  of  this  section  are  based  essentially  on  two  concepts: 
expansion  of  the  particle  trajectories  and  of  the  self- consistent  fields  in  the 
electron-ion  mass  ratio,        and  separation  of  the  Maxwellian  ion  distribution 
at   X  =  -  00    into  velocity  ranges  which  exhibit  entirely  different  trajectories 

2 

We  introduce  the  parameter   e     =  m    /m  More  precisely,    we  keep 

all  macroscopic  parameters  fixed     including  the  Alfven  speed  which  fixes    m    . 

2  2 

and  take   m     ~e    .      We  have  L     ~    1     L     ~  e    ,     and   L  ~  e.      An  important  point 

is  that  the  charge  separation  field  is  assumed  to  scale  as    l/e.     Because  of  this, 

the  motion  of  the  electrons  is  a  modification  of  the  conventional  guiding  center 

motion.      In  particular,    the  spiraling  frequency  is  not  the  gyro  frequency, 

2  1  /  2         -.- 

(J  =  eB/m    ,    but  is    [u     +  Lj(d/dx)(E    /B)]  =  u''',    the  adiabatic  invariant  is 

^2        .,  '  ^  ^2 

mv    /2u''~    rather  than  the  magnetic  moment    mv    /2cj.      The  guiding- center  drift 

velocities  m  the   x   and  y   directions  are  given  by 

m       dv 

v       =    -   E    /B  ,  u       =     E    /(B    +    -^    -— ^)     .  (27) 

-  X  -  y  e         dx 

To  lowest  order,    the  ions  see  only  the  large  electrostatic  field   E      (if 
they  are  in  the  body  of  the  shock     away  from   x  -  +  oo).      Ions  of  low  energy  could 
be  trapped  between  two  potential  peaks  of   E   .,      The  number  of  trapped  ions  in 
each  valley  must  be  known  in  order  to  compute  a  self- consistent  solution..     This 


Some  analysis  of  the  problem  on  this  basis  has  also  been  made  by  M.    Rosen- 
bluth  and  C     Longmire.     Private  communication. 
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is,    however,    misleading.  We  know,    from  the  exact  equations  (cf.    equa- 

tion (25)),    that  all  ions  must  have  come  from  the  Maxwellian  at   x  =  -  oo    and 
must  eventually  proceed  to  x  =  +  oo .      It  is  therefore  necessary  to  go  to  higher- 
order  approximations  in  the   e-expansion  (which  allow  the  ions  to  pass  through) 
to  determine  the  ion  states  which  are  inserted  m  the  lowest-order  solution. 

This  problem  takes  on  entirely  different  orders  of  complexity  depend- 
ing on  whether  all  the  ion  trajectories  are  monotone  in  x   or  whether  some  of 
the  orbits  loop.     In  the  second  case,    the  density  and  current  contribution  of  an 
ion  at  a  point  is  the  sum  of  the  individual  contributions  from  each  occasion 
when  the  ion  passes  that  point.     If  the  shock  is  not  too  strong,    then  an  ion 
started  at   x  =  -  oo    with  the  mean  fluid  x-velocity  and  no  y-velocity  (i.  e.    at  the 
peak  of  the  Maxwellian)  will  not  loop.     Neither  will  those  ions  loop  whose  ini- 
tial velocities  lie  within  a  certain  circle   C      m  velocity  space  surrounding  the 

2         2  2 

mean  ion  velocity,     (u  -  u   )     +  v      <  R    .      Consider  the  following  modification 

o  o 

of  the  problem  originally  posed.     There  are  assumed  to  be  no  ions  at   x  =  -  oo 

outside  the  circle   C    ,    the  distribution  is  a  Maxwellian  which  is  cut  off  at    R   . 

o  o 

To  lowest  order  in  the  expansion  parameter   e,    this  modified  problem  can  be 
solved  self- consistently.     The  result  turns  out  to  be  qualitatively  the  same  as 
that  described  m  Section  4.     There  are  periodic  solutions  and  pulses,    but  no 
shocks . 

An  ion  lying  outside  the  circle    C      may  be  reflected  back  by  the  elec- 
tric field   E    .      The  ions  can  be  classified  according  to  the  number  of  reflec- 

X 

tions  they  make  against  the  first  potential  barrier  which  they  meet  until  they 

pass  over  it  (as  eventually  they  must).     It  is  found,    to  lowest  order  in   e,    that 

after  having  overcome  the  first  barrier,    an  ion  will  not  be  turned  back  again 

for  at  least  the  distance   L/e.      Now  consider  a  concentric  circle    C      with 

radius    R,    >   R      within  which  the  ions  make  only  a  fixed  number  of  bounces  (the 
1  o  -^ 
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number  is  bounded,    independent  of  ^).^      This  problem  can  also  be  solved 
self- consistently  to  lowest  order  in   c.     The  solution  starts  as  a  pulse  and 
then  becomes  a  periodic  wave  (Fig.    7).     The  reason  for  the  distinction  be- 
tween the  first  wave  and  all  the  suc- 


B-B 


Fig.  7 


ceeding  waves  is  that  only  within  the 
first  are  any  ions  reflected  back. 

It  is  intuitively  appealing  to 
interpret  these  "collisions  of  ions  with 
»  the  electric  field"  as  being  similar  in 

effect  to  conventional  intermolecular 

collisions  in  an  ordinary  gas.     Consider 

2        2  2 

a  circle   (u  -  u   )    +v     »  R      on  which  some  ions  pass  through  the  potential  bar- 
o 

rier  without  reflection  while  others  pass  through  only  on  the  second  attempt. 
It  is  evident  that  this  circle  may  be  enormously  deformed  by  the  time  it  ulti- 
mately gets  entirely  past  the  potential  barrier.     After  passing  over  an  infinite 
oscillatory  wave  train,    it  would  not  be  surprising  if  this  circle  becomes  dense 
in  a  two-dimensional  domain.     A  smaller  circle,    on  which  there  are  no  reflec- 
tions at  all,    will  be  only  slightly  deformed  after  passing  through  a  wave.     This 
could  be  termed  a  "small-angle"  collision  of  the  ions  with  the  electromagnetic 
field  and  is  a  much  slower  process. 

It  must  be  noted  that  the  expansion  of  a  particle  trajectory  in  e    can- 
not be  expected  to  be  valid  uniformly  for  an  indefinite  length  of  the  trajectory. 
In  particular,    the  result  shown  in  Fig.    7  is  not  valid  for  distances  on  the  order 
of  1/e  wavelengths  from  the  first  potential  peak.     This  is  about  one  ion  phase 
length  (equation  (23))  and  is  the  distance  the  ions  would  go  before  being  affected 


To  lowest  order  in   e   the  circles    C      and  C      are  independent  of    e. 

o  1 
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appreciably  by  the  magnetic  field.  Consequently,  it  is  not  surprising  that  the 
second  type  of  damping  (termed  "small-angle"  above)  does  not  emerge  to  this 
order  m   e. 

An  interesting  group  of  ions  can  be  found  which  makes  on  the  order 
of   1/e  bounces  before  passing  through.     During  this  process,    they  gain  enor- 
mous energy.     Their  effect  is  hard  to  assess  quantitatively  since  each  such 
particle  contributes  a  large  damping  effect,    but  there  are  few  such  particles 
m  a  Maxwellian  distribution.     The  importance  of  even  a  very  few  such  highly 
energetic  ions  may  be  great  m  thermonuclear  applications.     This  mechanism 
is  reminiscent  of  the  Fermi  mechanism  for  accelerating  cosmic  particles, 
except  that  m  this  case  it  is  the  electric  field  which  is  dominant. 

It  is  clear  that  the  determination  of  the  width  of  a  shock  of  the  type 
envisioned  here  is  a  complex  matter.     The  result  depends  very  closely  on  the 
precise  shape  of  the  distribution  at   x  =  -  oo    (which  m  practice  could  be  rather 
arbitrary,    as  for  example  the  residue  left  by  an  earlier  shock!".     There  exists 
the  strong  possibility  that  the  greater  part  of  the  damping  occurs  across  the 
first  wave,    in  which  case  this  could  be  taken  as  a  practical  measure  of  the 
shock  width  even  though  the  final  decay  might  be  extremely  slow.     Such  a 
shock  would  be  observed  as  a  discontinuous  front. 

It  IS  possible  to  analyze  the  ultimate  decay  to  the  state  behind  the 
shock  by  utilizing  certain  stochastic  approximations       It  is  important  to  real- 
ize that  the  damping  mechanism  is  basically  nonlinear  and  cannot  be  expected 

to  produce  an  exponential  decay.     One  can  estimate  the  ultimate  decay  to  be  at 

1  /4 
a  rate    1 /x        ..      However,    if  one  lets   e   approach  zero  first  and  then  lets   x 

1/2 
approach  infinity,    a  decay  at  the  faster  rate    1 /x  is  plausible.     The  damp- 

ing length  scale  will  also  depend  on  the  strength  of  the  shock,    and  one  can 
estimate  it  as  becoming  exponentially  large  as  the  strength  approaches  zero. 
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For  shocks  which  are  stronger  than  a  certain  critical  value  beyond 
which  the  fluid  theory  of  Section  4  breaks  down,   there  are  indications  that, 
provided  they  exist,    there  are  at  least  two  stages  in  the  transition.     In  one 
stage  there  are  oscillations  whose  effective  wavelength  becomes  much  smaller 
by  a  factor   e"    which  brings  the  wavelength  below  the  gyro  radius  of  the  elec- 
tron     and  charge  separation  becomes  an  important  feature  of  the  problem. 
The  roles  of  electrons  and  ions  with  regard  to  heating  are  reversed  in  this 
stage.     It  is,    however,    impossible  to  estimate  qualitatively  whether  or  not  the 
ions  (heated  m  the  other  stage)  end  up  hotter  than  the  electrons. 


7.  A  RELATED  TRANSIENT  PROBLEM 

The  question  of  the  approach  to  equilibrium  behind  a  shock  is,    m 
principle,    better  answered  by  studying  the  transient  problem  (e.  g.    of  a  piston 
moving  into  the  plasma)  than  by  studying  the  stationary  problem.     On  the  other 
hand,    the  transient  problem  is  more  difficult.     In  this  section  we  discuss  the 
solution  which,    roughly  speaking,    describes  the  compression  wave  produced  by 

a  slight  piston  motion.     For  simplicity  we  take  the  case  of  equal  ion  and  elec- 

2 
tron  masses,    thereby  eliminating  the  charge  separation  field,     E    .  "     The  pres- 
sure is  taken  to  be  zero;     consequently  we  are  free  to  consider  these  equations 
as  representing  a  fluid  or  a  self- consistent  particle  stream.     The  x-component 


Arguments  in  which  the  electron  gyro  radius  dominates  are  detailed  in  Colgate, 
S.  ,    A  Description  of  a  Shock  Wave  m  Free  Particle  Hydrodynamics  with  Inter- 
nal Magnetic  Fields.    UCRL-4829,    (1957), 

2 
Time-independent  problems  of  this  type  have  been  treated  by  M.    Mittelinan 

and  M.    D.    Kruskal.     Private  communication.     Also  see  paper  by  Adlam  and 

Allen  in  these  Proceedings. 
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of  velocity  is  the  same  for  ions  as  for  electrons;    the  y- components  are  nega- 
tives,    V     =  -  V   .      The  number  of  equations  is  thereby  reduced  to  that  in  a 
single  fluid  theory.     For  linear,    small  amplitude  departures  from  a  state  of 

rest  (as  in  Fig.    1,    but  with  the  unperturbed  state   u     =0   and   E       =0   and 

o  yo 

n  =  n   ),    the  equations  are 
o 

at        9x 

m-—    +    eB    V    a    0 
dt  o 

m-^    -     e(E     -  uB    )     =0  (28) 

9t  y  o 

9B 

-r —    +    uen   V      =0 

9x  o 

9E 

9x  9t 

It  is  possible  to  reduce  this  system  to  a  single  equation  for  E    . 

2  2  4 

9    E  ^9    E  ,    9    E 

-f    -    A'-f    =    L^-V^      ■  (29) 

9t  9x  9x   9t 

A  is  the  AlfvSn  speed  (13)  and  L  is  essentially  the  same  reference  length  as  in 
(22),  taking  note  that  m  =  m  and  that  the  shock  speed,  in  so  far  as  it  makes 
sense  to  speak  of  it  in  this  problem,    would  be  A. 

It  is  convenient  to  impose  as  a  boundary  condition  at   x  =  0   the  sudden 
application  of  a  constant  transverse  electric  field   E        rather  than  a  piston 
condition.     If  one  imagines  an  originally  quiescent  cylindrical  plasma,    then  a 
suddenly  applied  tangential  electric  field  corresponds  to  a  conventional  pinch. 
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Equation  (29)  with  the  given  initial  and  boundary  conditions  can  be 
solved  by  a  Laplace  transform  with  respect  to  t.     The  inverse  of  the  trans- 
form cannot  be  found  in  closed  form,    but  its  asymptotic  behavior  for  large  t 
can  be  determined  by  the  method  of  steepest  descents. 


E   (x, t)    ~    E 

y  yo 


00 


(s)ds 


1/3 

(-)        a 
^3 


where 


1         x-At 


and    A.{s)     is  the  Airy  function,    defined  by 

oo 

A.(s)     »    -         \     cos(-y     -  sy)  dy 
1  '^         J  3 


This  result  represents  a  wave  train  moving  at  the  Alfv^n  speed. 
Figure  (8)  shows  the  qualitative  behavior  of  E     plotted  against  a,    for 


tEy(X.t) 


yo 


1. 


Fig.  8 


The  presence  of  the  cube  root 
of  time  in  the  denominator  of  the  scale 
parameter  a   shows  that  the  width  of 
the  wave  is  increasing  as  the  cube  root 
of  time.     The  height  of  the  wave  drops 
off  exponentially  ahead  of  the  front, 
C3f  ■  0   or   X  *  At.      For  negative  a,    the 


curve  oscillates,    and  approaches  unity.    The  amplitude  of  the  oscillations  de- 

-3/4 
creases  like  a 
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This  result  is  not  easily  compared  with  the  previous  ones       By  anal- 
ogy with  ordinary  gas  dynamics,    one  would  guess  that  the  spreading  of  the 
wavefront  is  due  to  the  linearization       However,    we  know  the  solution  of  the 
appropriate  non-linear  equations  in  the  steady  state,    and  there  is  no  corre- 
sponding steady  shock.     Also,    the  damping  mechanism  is  itself  nonlinear,    so 
the  damped  shape  of  this  linearized  wave  may  not  be  representative  of  the 
answer  to  the  preceding  problem.     In  this  context,    we  can  refer  again  to  the 
gas  dynamical  shock.     The  solution  to  a  free-flow  initial  value  problem  (i.  e. 
with  no  collisions)  widens  indefinitely  even  when  of  finite  amplitude.     Colli- 
sions are  required  to  validate  the  fluid  equations  (which  imply  steepening)  as 
ell  as  to  supply  dissipation  (the  conventional  widening  mechanism). 


w 


8       CONCLUSIONS 

The  object  of  study  has  been  a  medium  strength  shock  (measured  by 
M'''    or  by  the  magnetic  field  ratio).     The  conclusion  that  the  major  portion  of 
the  heating  is  applied  to  the  ions  follows  immediately  from  the  single  assump- 
tion that  the  electrons  are  adiabatic.     The  goal  of  this  entire  complex  self- 
consistent  theory  is  only  to  verify  this  assumption  as  well  as  the  implicit  one 
that  the  shock  transition  zone  is  small  enough  to  fit  mside  the  laboratory. 
The  most  interesting  case,    producing  a  very  large  temperature  rise,    is  when 
the  plasma  energy  m  front  of  the  shock  is  a  small  fraction  of  the  inagnetic 
energy  (/3     «  1)       However,    this  is  exactly  the  case  when  the  wave  tram  is 
largest.     If  quantitative  computations  verify  that  a  significant  part  of  the  en- 
tropy increase  takes  place  across  the  first  wave,    then  the  macroscopically 
small  length   L   (equation  (22))  can  be  taken  as  the  effective  width  of  the  shock. 
If  this  is  the  case,    it  should  be  possible  to  obtain  thermonuclear  temperatures 
from  a  relatively  cold  plasma  by  the  passage  of  a  single  shock  wave.      Whether 
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experiment  or  theory  will  be  the  first  to  definitely  verify  or  refute  this  possi- 
bility remains  to  be  seen. 

REFERENCES 

Many  of  the  results  described  in  this  paper  were  originally  presented 
at  Sherwood  meetings;    e   g     in  Berkeley,    February.    1957,    TID-7536  (Pt.    1), 
pp.    32-38;     in  Washington.    February,    1958,    TID-7558     pp.    307-310,    311-312. 
Also  see  Morawetz,    C.    S.  ,    Magneto- hydrodynamic  Shock  Structure  Using 
Friction,    NYO-8677,    New  York  University,    Institute  of  Mathematical  Sciences, 
(January  1959). 


-  28 


PROBLEMS  IN  MAGNETOHYDRODYNAMIC  STABILITY 


J.    Berkowitz,    H.    Grad,    H.    Rubin 


1 


vacuum 


1.     INTRODUCTION 

Our  primary  concern  is  the  hydromagnetic  stability  of  free-boundary 

2 
equilibria.         These  are  configurations  in  which  a  perfectly  conducting  field- 
free  fluid  is  separated  at  an  interface  (surface  charge  and  current)  from  a  vac- 
uum electromagnetic  field   (E,  B),    (Fig.    1).     Two  examples  were  treated  by 

3 
Kruskal  and  Schwarzschild:       one,    a  conducting 

fluid  supported  against  gravity  at  a  plane  inter- 
face by  a  uniform  magnetic  field  (Taylor  insta- 
bility);   the  other,   an  infinite  cylindrical  con- 
ducting fluid  surrounded  by  an  azimuthal  mag- 
netic field  (pinch  instability).     In  full  generali- 
ty,   we  consider  a  free-boundary  equilibrium 
with  an  interface  of  arbitrary  shape  in  the  pres- 
ence of  an  arbitrary  conservative  force  field.     The  configuration  of  external 
fixed  conductors  is  arbitrary.     The  external  magnetic  constraints  may  be 
either  constant  fluxes  or  constant  currents.     The  conducting  fluid  may  be  taken 
incompressible  or  compressible  with  any  equation  of  state,    dissipative  (viscous 


Grateful  acknowledgment  is  made  to  A.   A.    Blank,    K.    O.    Friedrichs,    and  H. 
Kranzer  for  their  assistance. 
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or  heat  conducting),    non-dissipative,    or  even  governed  by  the  Boltzmann  equa- 
tion,   with  or  without  collisions.     Stability  is  considered  both  with  respect  to 
small  amplitude  and  finite  amplitude  perturbations. 

It  is  possible  to  solve  the  problem  in  such  general  terms  because  the 
stability  problem  separates  into  an  electromagnetic  problem  which  can  be 
treated  as  essentially  static  if  displacement  current  is  ignored,    and  a  fluid 
problem  for  which  stability  is  treated  by  conventional  thermodynamic  consider- 
ations.     The  stability  criterion  is  expressed  m  terms  of  geometric  properties 
of  the  interface  alone  and  is  independent  of  the  constraints  or  of  the  equations 
which  govern  the  fluid.     A  number  of  absolutely  stable  equilibria  are  given 
(Section  6), 

Part  of  this  work  has  been  extended  to  the  case  of  conducting  fluid 

and  magnetic  field  which  are  intermingled  instead  of  separated  at  an  inter- 

2 
face.         Certain  aspects  of  this  more  general  problem  are  treated  here  (Sec- 
tion 10),    and  a  class  of  stable  axially  symmetric  configurations  is  described. 
The  more  difficult  problem  of  stability  in  the  presence  of  periodically  oscillat- 
ing magnetic  fields  is  solved  for  a  special  geometry  (Section  11). 

For  a  Hamiltonian  system  with  finitely  many  degrees  of  freedom,    the 
type  of  argument  used  in  this  stability  analysis  yields  categorical  answers. 
For  such  a  system,    an  absolute  minimum  of  a  monotone  potential  represents 
a  stable  equilibrium  in  the  sense  that  a  small  perturbation  will  never  cause 


The  methods  and  results  described  in  this  paper  were  originally  presented 
at  Sherwood  meetings  m  Princeton,    October  1954  (WASH- 184,    p.    144),    in 
Berkeley,    February  1955  (WASH- 289,    p.    115),    in  Princeton,    October  1955 
(TID-7503,    p.    238,    p.    241,    and  p.    247),    and  m  Gatlinburg,    June  1956  (TID- 
7520,    p.    373). 
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Bernstein,    I.   B.  ,    E.   A.    Frieman,    M.    D.    Kruskal  and  R.    M.    Kulsrud,    An 

Energy  Principle  for  Hydromagnetic  Stability  Problems,    NYO-7315,    Matter- 
horn  Project,    Princeton  University,    Princeton,    N.    J. 
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more  than  a  small  deviation  from  equilibrium,    and  a  finite  perturbation,    a 
correspondingly  finite  deviation.     If  the  equilibrium  state  corresponds  to  a 
relative  minimum  of  the  potential,   there  exists  a  maximum  perturbation  of 
the  energy  below  which  the  system  remains  in  some  neighborhood  of  the  equili- 
brium.    The  existence  of  the  relative  minimum  (but  no  estimate  of  the  maxi- 
mum stable  perturbation)  can  be  ascertained  by  examining  the  second  varia- 
tion of  the  potential  function.     From  the  second  variation,    we  may  obtain 
more  than  a  stability  criterion:    the  eigenvalues  computed  with  the  kinetic 
energy  as  unit  form  are  the  frequencies  (stable)  or  the  growth  constants  (un- 
stable) of  small  amplitude  disturbances. 

For  a  fluid,    the  infinity  of  degrees 
of  freedom  introduces  complexities.     Con- 
sider the  example  of  a  container  of  water 
under  gravity  (Fig.    2).     There  exists  a  po- 
Fig.  2  tential  and  an  energy  integral   K+P   (kinetic 

plus  potential).     In  Fig.    2a,    P     is  the  abso- 
lute minimum  of  P.      For  an  energy  perturbation  e,    we  have   K  +  P  »  P     +  e. 
We  conclude  without  further  reference  to  the  equations  of  motion  that  for  all 
time 


0  <  K  <  e 


P    <  P  <  P     +  e 
o  —      —     o 


(1) 


The  first  inequality  states  that  the  kinetic  energy  is  never  large.     The  second 
states  that,    in  some  sense,    the  configuration  never  deviates  very  much  from 
equilibrium;    for  example,   that  no  more  than  a  certain  small  amount  of  water 
can  exceed  a  given  height.     It  is  intuitive  that  the  same  conclusion  applies  to 
Fig.    2b,    viz.   that  for  a  given  e   only  a  certain  amount  of  water  can  be  lost  over 
the  rim.     This  conclusion  does  not  follow  from  energy  considerations  alone  and 
the  inequalities  (1)  cannot  be  rigorously  proved  even  for  restricted  values  of  e. 
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In  principle,    it  may  be  possible  for  a  fluid  particle  to  gain  large  kinetic  energy 
after  passing  over  the  rim  and  pass  this  energy  back  to  the  water  in  the  con- 
tainer.    Because  of  this  almost  insurmountable  mathematical  difficulty,    the 
stability  of  a  fluid  system  with  respect  to  finite  amplitude  perturbations  is 
treated  in  two  parts,    the  first,    mathematically  precise,    will  show  that  a  given 
equilibrium  minimizes  an  appropriate  potential  among  all  neighboring  states, 
the  second  will  be  a  heuristic  evaluation  of  this  result  m  the  light  of  some  "prac- 
tical" criterion  for  stability. 

The  point  of  this  variational  approach  is  to  eliminate  almost  all  refer- 
ence to  the  equations  of  motion  or  their  solutions,     We  use  only  the  integrals 
and  constants  of  the  motion,    energy  and  others       Although  we  would  like  to  com- 
pare the  potential  of  the  equilibrium  state  with  the  potentials  attainable  by  possi- 
ble motions  of  the  system,    we  compare  it  with  the  potentials  of  all  mathemati- 
cal configurations  compatible  with  the  chosen  constraints  (i.  e.  ,    the  known  con- 
stants of  the  motion).     The  manifold  of  stable  solutions  could  be  incomplete  if 
there  exist  constants  of  the  motion  which  are  not  recognized;    also,   the  method 
may  be  too  conservative  if  neighboring  configurations  exist  which  are  not  dy- 
namically accessible  as  solutions  of  the  equations  of  motion. 

A  positive  second  variation  of  the  potential  for  all  admissible  infini- 
tesimal perturbations  is  evidence  for  stability.     This  finding  makes  plausible 
the  existence  of  a  relative  minimum  but  does  not  insure  it  and  is,   therefore, 
a  less  conservative  criterion. 

An  entirely  different  kind  of  stability  analysis  is  to  linearize  the  equa- 
tions of  motion  and  separate  an  exponential  time  factor  (normal-mode  analysis). 
It  is  often  possible  to  relate  the  second  variation  formulation  of  the  stability 
problem  to  the  normal  mode  analysis  as  m  the  case  of  finitely  many  degrees 
of  freedom.     The  existence  of  an  unstable  mode  is  better  evidence  of  instabil- 
ity than  the  mere  existence  of  a  perturbation  which  makes  the  second  variation 
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negative,   but  the  normal  mode  analysis  is  no  more  trustworthy  in  predicting 
stability  than  the  second  variation. 

The  stability  of  fluid  magnetic  equilibria  will  be  considered  from  all 
three  viewpoints-    minimization  of  an  appropriate  potential  with  regard  to 
finite  perturbations,    analysis  of  the  second  variation,    and  normal  mode  analy- 
sis.    For  brevity,    we  touch  only  upon  the  major  points,     A  detailed  presenta- 
tion will  be  given  elsewhere. 

2.     CONSTANTS  OF  THE  MOTION 

For  the  electromagnetic  field  we  use  Maxwell's  equations  (in  ration- 
alized M.  K.  S.    units)  without  displacement  current. 


^    +    curl  E    =    0    ,  S  (2) 


At  the  interface,    moving  with  fluid  velocity  u,    we  have  the  boundary  condi- 
tions 

B       =    0     ,  (E  +  uxB)^     =    0     ;  (3) 

n  t 

the  subscripts   n,  t   denote  normal  and  tangential  components  respectively. 
Only  u     actually  enters  in  (3).     At  a  fixed  conducting  boundary,    we  have 

B       =    constant  in  time,  E      =    0.  (4) 

n  t 

The  electromagnetic  force  JxB   is  normal  to  the  interface  and  exerts 

2  , 
a  pressure  of  magnitude   B    /2/j.      The  electric  force   qE  and  also  the  dis- 
placement current  are  relativistically  small;    moreover  the  two  must  be 
dropped  or  kept  together  to  keep  the  electromagnetic  momentum  equation  in 
conservation  form.     The  system  (2)  and  the  boundary  conditions  (3),    (4)  are 
Galilean  invariant  and  are  even  preferred  to  the  Lorentz  invariant  Maxwell 
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system  for  use  here  with  non-relativistic  fluid  equations. 

Under  the  boundary  conditions  (3),    (4/  the  solution  of  the  system  (2) 
is  determined  at  each  instant  by  the  position  of  the  interface  without  regard 
to  its  velocity. 

If  B     =  0   at  all  boundaries,    fixed  or  fluid,    the  magnetic  field  is  de- 
termined by  a  finite  number  of  parameters,    either  fluxes  or  m.  m.  f.  '  s  (cur- 


rents) 


r 

J      =         I    B'  dS    ,       or       I       =    -  B-  dx      .  (5) 

-^r  J  r         M       J 

S  C 


The  magnetic  energy  is  then  given  by 

M    =       \^  B^dV    =    /      -  I   J 

9.11  ^ '    9       r^-J- 


2\x  ^ — '  2    r-r 

r 


(6) 


=   ^\l    L       I       =    'Y\1    X       $      . 

^ '  2    r     rs    s         ^ — '  2  -^r    rs  -^s 

r,  s  r,  s 

where  the  inductance  matrix  L        and  the  susceptance  X        are  defined  by 

rs  rs 


^=>LI,            I=/X(|). 
_LjP         ^. — ,     rs   s  r         ^ rs  -^s 

From  the  integrated  form  of  the  first  equation  of  (2),    we  have 


Blank,    A.    A.,    H.    Grad  and  K.    O.    Friedrichs,    Theory  of  Maxwell's  Equa- 
tions without  Displacement  Current,    NYO-6486,    New  York  Umversity,    Insti- 
tute of  Mathematical  Sciences  (1957), 


-  34 


r  d 

— i    =   -f-  BdS 

dt  dt 


(E  +  uxB)  •  dx    =    0      ,  (8) 


S   (t)  C   (t) 

r  r 

where  C      is  the  boundary  of  S   .      We  conclude  that  the  fluxes  are  constants 

r  '^  r 

of  the  motion.     An  example  is  a  toroidal  plasma  surrounded  by  a  vacuum  mag- 
netic field  contained  within  a  concentric  copper  torus.     Two  fluxes  {or 
m.  m.  f.  '  s)  completely  determine  the  instantaneous  magnetic  field.     In  any  de- 
formation of  the  plasma,    the  two  fluxes  are  constant. 

More  generally,    when  magnetic  lines  intercept  a  part  of  the  boundary 

(B     ^  0),    the  magnetic  field  is  determined  by  the  given  stationary  values  of 

B      as  well  as  certain  fluxes   (j)    ,     these  fluxes  again  are  constants  of  the  mo- 
n  -*-r  ^ 

tion. 

Constant  current  generators  may  be  introduced  by  violating  the  condi- 
tions (4).     For  example,    we  may  introduce  slits  in  the  fixed  conductors  and 
apply  suitable  potentials  across  them.     With   E    =  0    elsewhere,   the  potential 
across  the  slit  will  be  equal  to  the  line  integral  in  (8/,    destroying  9  =  constant 
but  permitting  I     =  constant.     More  generally,    we  could  distribute  the  voltage 
arbitrarily  around  the  conductor  so  long  as  the  total  e,  m.  f.   is  sufficient  to 
keep  the  current  constant.     Such  a  distributed  e.  m.  f,    can  be  consistent  with 
the  boundary  condition  B     =  constant   and,    as  before,    the  magnetic  field  is 

uniquely  determined  by  the  given  values  of  B      and  the  fixed  currents    I      (in- 

n  r 

stead  of  d)    ).         This  result  holds  even  for  a  mixture  of  constant  current  and 
-^r 

constant  flux  constraints.     In  most  of  the  subsequent  analysis,    we  shall  take 
the  general  case,    there  exist  n   independent  fixed  circuits,   the  first   m   of 
constant  current  type  and  the  remainder  of  constant  flux  type; 

See  footnote,    p.    34. 
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I     =  constant,     (r  *  1,  2,  ...  ,  m),        q)     =  constant,     (r  =m+l,  .  .  .  ,  n) 


(9] 


The  rate  of  change  of  magnetic  energy  in  a  moving  domain  V(t)    as 
computed  from  (2),    (3),    (4)  is       :'      '   ' 


dM    _     d 
dt      ~    dt 


f  B^dV 


^  B^u  -  -  ExB 


dS 


V(t:f 


S(t) 


i)     ^  B   u-dS    -       ()     -  ExB-dS    , 
J      2a(  J     m 

S  (t;  S  (t)  ■ 


(10) 


where   S     denotes  the  interface  and   S     the  fixed  boundary;    the  Poynting  vec- 
tor  contribution  on   S     reverses  the  sign  of  this  integral.     It  can  be  shown  that 
only  the  constant  current  circuits  of  (9)  contribute  to  the  second  integral; 

namely,       ,     ....-,      ...........       ^ .  ^  .    ■.    .   ■ 


m 


m 


^     \   ExB-    dS    =    -X^I    ^    =    -4Z^I     ^ 
u  ^—-'    V   dt  dt  ^—^    r  -r 

J  r=l  r=l 


(11) 


If  we  set 


M'''     =    M 


Sz:  :  .  V  :■ 


m 


Z'j 


r  =  l 


r  -^r 


(12) 


we  have  in  the  general  case,    (9) 


dM'" 
dt 


^  B^u-  dS 
2m 


(13) 


i: 
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In  the  special  case  B     =  0   on  all  boundaries  we  have  from  (6) 


n 


n 


m 


M 


*     _ 


Z  |i  i  -Z 


^.   2  >^r         ^—^   2    r  ^r  ■ 
r=m+l  r  =  l 


(14) 


If  the  fluid  IS  incompressible,   we  have  the  equations 


div  u    =    0  ;       1^   +   (u-Viu   +-Vp   +  V0     =    0 
9t  p 


(15] 


from  which  we  obtain  the  rate  of  change  of  energy, 


dt      J     2  ^"   ^^    -    - 
V(t)  S(t) 


I     (p  +  p  0)u' 


dS 


(16) 


where  p   is  the  mass  density,    u  is  the  fluid  velocity,    p  is  the  pressure,   and 
0   IS  an  external  potential  energy  per  unit  mass.     For  any  other  kind  of  fluid, 
taking  heat  flow  to  be  zero  across  the  interface  (radiation  is  neglected),    we 
have 


_d 

dt 


pie  +1  u^  +  0)dV    =    - 


u.P.dS. 
1    ij      J 


(17) 


V(t) 


S^(t) 


where   e  is  the  internal  energy  per  unit  mass,   and   P. .   is  the  stress  tensor. 
The  electromagnetic  stress  is  normal  to  the  interface  and  must  be  counter- 
balanced by  the  fluid  stress  (otherwise  the  acceleration  would  be  infinite); 
hence,    the  fluid  stress  reduces  to  its  normal  component    p     and 


p^    =    B    /2tx  (on  S^)    , 


(18) 


whence 


dK^dt* 
dt       dt 


p  u-dS 

n 


(19) 
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where 


K 


1  2 

-  pu   dV 

2  ^ 


(20) 


J 


d* 


p(e  +  0)dV 


(21) 


From  (13),    observing  that  the  normal  to  S      is  oppositely  oriented  from  the 
magnetic  and  fluid  sides     we  obtain  the  energy  integral 


K  +  6*  +  M"     =    constant. 


(22) 


In  the  incompressible  case  we  set    e  =^  0    in  <£'"    and  have  an  additional  integral, 
the  total  volume  of  the  fluid,  \ 


J 


dV     =    V 


constant. 


(23) 


3,     EQUILIBRIUM.     SEPARATION  OF  THE  POTENTIALS 

In  a  non-dissipative  fluid     tnechanical  equilibrium  is  characterized 
by  the  pressure  balance, 


1 


Vp  +  V0   =  0 


(24) 


This  implies  that    p  and  p    (and  hence  all  other  thermodynamic  quantities*  are 
constant  on  the  surfaces  of  constant    0        If  the  fluid  is  in  thermodynamic  equili- 
brium,   the  temperature  is  constant  throughout.     Employing  the  Gibbs  free 
energy  per  unit  mass,     F  =  e  -  Trj  +  p/p.     we  obtain    dF  =  -  r/dT  +  —  dp   (r)  is  the 
entropy  per  unit  mass),    and  from  (24),    with   T    constant. 


F  +  0     =    constant. 


(25' 
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For  an  incompressible  fluid,    we  have 

p  +  p0     =    constant,  (26) 

In  equilibrium  (scalar  pressure),    the  boundary  condition  (18)  takes 
the  form 

p    =    B^/2m  (27) 

For  an  arbitrarily  given  interface  there  exists  a  unique  equilibrium  fluid  con- 
figuration on  one  side  and  a  unique  magnetic  field  on  the  other.     The  boundary 

condition  (27)  singles  out  a  special  class  of  equilibrium  shapes.     In  the  sim- 

,  2 

plest  case,    0=0,     we  have   p  =  constant   in  the  fluid  and  B     =  constant  at  the 

interface.     Several  specific  examples  will  be  given  later.     The  problem  of  de- 
termining such  equilibria  is  identical  to  the  free-boundary  problem  for  irrota- 
tional  incompressible  fluid  flow   (B    corresponds  to  fluid  velocity)  where  the 
condition  is  that  the  absolute  velocity  is  constant  at  the  unknown  surface. 

We  now  define  the  potentials    P  in  a  way  which  separates  magnetic 
and  fluid  effects.     In  a  specific  equilibrium  state,    the  fluid  pressure  is  a  func- 
tion of  position   X    defined  by  the  implicit  relation  (25)      F(p(x' ,  T)  +  0(x)  =  con- 
stant.     We  suppose   0(x^   to  be  defined  in  all  space.     The  definition  of   p(x)    may 
then  be  extended  into  the  vacuum  domain  using  (25);    the  extended  function  is 
called   p   (x).     Referring  to  (21)  we  put 


o 

r 


f  ^ 


p   dV     =      '      (pe  +  p0  +  p    ;dV  (28) 


o  I        ■  ■  *  o 


J        ^  J 

V  V 

f  f 


P_     =    M""  +     /       p^dV     =      /       (-;7::B"    +   P„'dV  -/      I    \       (29) 


m 
1    „2 


m  Jo  2[x  o '  r— r 

m 

where  V„  and   V       denote  the  fluid  and  vacuum  domains  respectively.     The 
f  m  i-  J 
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integral  of  p      has  been  added  to  both  the  fluid  and  vacuum  potentials.     Since 

V.  +  V       IS  constant  during  the  motion,    we  have  merely  changed  the  value  of 
f  m  s,  .  J  & 

the  constant  in 

K  +  P.  +  P        a    constant.  (30) 

f         m 


4.     VARIATION  OF  THE  FLUID  POTENTIAL 

The  stability  of  the  fluid  is  treated  by  purely  thermodyn.amic  argu- 
ments.    We  show  for  a  fluid  satisfying  any  conventional  equation,    dissipative 
or  not,    that  motion  proceeding  from  an  initial  thermodynamic  equilibrium 
state  can  only  increase  the  potential    P         Also,    for  a  non- dissipative  fluid 
with  no  external  potential   0,    a  motion  proceeding  from  a  state  of  dynamic 
equilibrium   (p   is  constant  but   T   is  not  necessarily  constant^;,    can  only  in- 
crease   P  . 

We  recall  that  the  total  entropy  of  a  system   S{(5,  V,  m^    is  a  homogene- 
ous first  order  function  of  the  three  arguments,    total  internal  energy,    volumej 
and  mass.     We  denote  the  conventional  densities  (per  unit  mass^  of  entropy  and 
internal  energy  by  -q,    e,     respectively,    and  introduce  the  densities  per  unit 
volume  of  entropy  ri  "^  p-q  si  nJT,    and  energy   e  =  pe  =  e/r^    where  t  ~  l/p   is 
the  specific  volume.     The  second  law  of  thermodynamics  is  equivalent  to  the 
statement  that  the  functions    e{ri,p':    and   e()7,  t)    are  convex  functions  of  their 
arguments  and  satisfy 

Tdr]    =   de  +  pdr        ;  Tdrj    =    de  -  Fdp  (31) 

At  a  given  point   x,     consider  the  function 

a(^,p)     =    e(^,p)    +    0p    -    Tji     ,  (32) 
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where  T     is  the  equilibrium  value  of  T.      From 
o 

9^    _    9e  _ 

— 7\      ~     — 7\       "1  -       i        -        1 

dr]  dr]  o  o 


(33) 


9«  9e  .  ^  . 

we  see  that  a  is  stationary  atT=T,     F  +  0=O.      The  function  a   is  convex 

o 

in  n   and  p    since  it  differs  from    e  by  a  linear  function.     We  conclude  that  a 
has  a  unique  minimum  at  values  of  rj   and  p    corresponding  to  T  =  T      and 
F  +  0  =0.      The  minimum  property  is  not  altered  by  adding  the  constant   p 
(at  the  given  point   x)   to  a,    and  we  have 

e  +  0p-Tn  +  p      >e+0p-Tr)+p       =p(F+0;=O.  (34) 

oo  o  o         ooo  oo 

Integrating  over  the  fluid  volume  we  obtain 
(^''  -  T   S  + 


°        J 


pdV>0=<5'-TS+         pdV  (35) 

o  o         o   o      i        o 


where    c'    is  defined  by  (20).     Since  the  total  entropy  S   cannot  decrease  in  any 
motion,    S  >  S   ,    we  conclude  that  motion  proceeding  from  an  initial  equili- 
brium can  only  increase  the  fluid  potential   P  . 

In  the  special  case  0=0   it  is  possible  to  relax  the  condition 
T  =  constant   of  thermodynamic  equilibrium^    and  by  integrating  with  respect 
to  mass  in  the  inequality 

e  +  pT-Tr]     >e     +pT     -Tr?       =F  (36) 

o  o      —      ooo         oo  o 

we  obtain  the  same  result  for  adiabatic  variations  from  a  dynamical  equili- 
brium (constant  pressure). 
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If  the  fluid  equations  of  motion  are  replaced  by  the  Boltzmann  equa- 
tion,   we  use  the  Boltzmann   H-function  instead  of  S, 


H 


f  logf  d?. 


(37) 


J 


where  the  molecular  distribution  function   f(?jX,  t)   is  a  function  of  velocity  C 
as  well  as   x  and  t.      In  equilibrium,     f  is  Maxwellian  and 


H 


1  ^ 


(38) 


where   R   is  the  gas  constant.     If  f  is  not  Maxwellian,    we  have 

H   >    -  ^  n  (39) 

where  r\(e,  p)   is  the  usual  thermodynamic  function  of  the  energy  and  mass 
density  associated  with   f. 


fd?     , 


e    = 


^C'fd? 


(40) 


From  the  Boltzmann  equation,    with  appropriate  boundary  conditions  at  the 
interface^    one  can  obtain, 


H(x,  t)dV    < 


H(x,  0)dV 


J 


For  an  initial  state  of  equilibrium,    we  have  from  (38),    (39), 


r7(x,  t)dV    > 


J 


RH(x,  t)dV    > 


RH(x,0)dV    =     1   n    dV  (41) 

o 


From  the  Boltzmann  equation  we  obtain  the  same  constraint  on   S  which  was 


42 


used  in  proving  the  previous  inequality.     Only  the  primitive  variables   p   and   e 
are  obtained  directly  from   f;     all  other  thermodynamic  quantities  are  defined 
as  the  conventional  functions  of  p    and   e.      The  fact  that  these  definitions  are 
arbitrary  and  generally  without  significance  for  large  deviations  from  equili- 
brium does  not  affect  the  validity  of  the  final  inequality  (35)  containing  the  un- 
ambiguous quantities   d''"    and  V.      The  parameter   p     is  also  unambiguous  as 

o 

the  pressure  in  a  state  of  equilibrium. 

We  have  shown  that  the  fluid  potential,     P  ,    is  an  absolute  minimum 
at  equilibrium.     For  completeness,    we  give  the  second  variation  of   P..      If  the 
variation  is  adiabatic,    computation  yields 


^\ 


[div(6x)]^pc^dV    ,  (42) 


2 

where   c     =  9p(r7,  p)/9p   and  the  result  is  formally  independent  of  0.     The  sec- 
ond variation  vanishes  for  any  incompressible  variation   div(6x)  =  0.      When 
0=0,    it  is  easy  to  see  that    P     is  unaltered  by  incompressible  variations  even 
if  they  are  finite. 

5.     VARIATION  OF  THE  MAGNETIC  POTENTIAL^ 

In  variational  analysis  it  is  conventional  to  define  the  variation    6B    of 
an  independent  function   B(x)   as 


6B     »    |-B(x,  e) 


e=0 


where  B(x)    is  imbedded  in  the  family  of  functions   B(x,  e)    so  that   B(x)  =  B(x,  0). 
We  extend  this  notation  m  not  insisting  that   €  =  0   in  (43).     For  the  variation 
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6x   of  the  interface,    we  set 

6x    =   |-x(S,  e)    , 

where  S   denotes  a  pair  of  coordinates  fixed  to  the  surface  as  it  varies. 

Since   div  6B  =  0,    a  vector  potential  exists  with   6B  =  curl  A.      The 

condition  B     =0   further  restricts    6B,    but  the  restriction  on   6B   (which  is 
n 

complicated  if  the  boundary  is  varied)  is  most  easily  written  in  terms  of  the 
potential-    the  condition   (A  +  Bx  6x)    =  0,    for  example,    will  guarantee  that 
B     =  0.      In  the  light  of  these  remarks  we  see  that  the  pre- Maxwell  equations 
(2)  may  be  interpreted  as  variational  equations  for  a  harmonic  (i.  e.  ,    irrota- 
tional  and  incompressible)  vector  field  where  t   correspcnds  to   e,    B  =  9B/9t 
to   6B,    u  to   6x,     E  to   -A,    etc.     We  shall  find  it  convenient  to  use  time  de- 
rivatives interchangeably  with  "variations".     In  particular,    we  shall  use  the 
"variational"  formulas 

r 

^      1     X(x,t)dV    =       I     ^dV    +        I     Xu-dS  (44) 


X(x,t)dV    =      J      l^dV    +  Xu- 


V(t)  V(t)  S(t) 


■^      I     A(x,t)-dS    =  |-|^  +  u  divA  j  -dS    +      (p    A*u-dx, 

S(t)  S(t)  C(t) 


(45) 


in  (44)    S(t)   is  the  boundary  of    V(t)  ,    and  in  (45)   C(t)   is  the  boundary  of  S(t). 
First  let  us  consider  the  magnetic  energy,     M  of  (6)  as  a  domain 


functional.     Applying  (44),    we  have 


6M 

■J 

-(B-B)dV    + 
M                             J 

V 

S 

^B^u    dS  (46) 

2/j 


m 


-  44 


The  relevant  potential  is    M'''    defined  by  (12),    for  which  we  obtain  after  some 
manipulation  the  first  variation 


6M'' 


(])    ^  B"u-dS 
2ju 


(47] 


The  same  formula  is  obtained  for  any  variation  of  currents,  fluxes,  and  of  the 
interface  subject  to  the  constraints  (9).  To  complete  the  first  variation  of  the 
vacuum  potential  (28),    we  obtain 


6  p   dV    =      d)      p  u  •  dS 


(48) 


V 


m 


In  equilibrium,    B    /2ju  =  p      on  the  interface,    cf.    (27),    and  on  adding  (47),    (48) 

we  conclude  that    6P       =0   and  the  potential  is  stationary      To  calculate  the 

m  -^ 

second  variation,    we  apply  (45)  to  obtain 


S^P 


m 


B-B  +  (u-Vh|b^)  I  u-dS 


1        2    . 
-r—  B    (u  +  u  divu)/  •  dS 

2/Li 


j    jp^u  +  u 


div  p  u  (  •  dS 


(49) 


Under  the  constraint  conditions  (9)  we  can  show  that 


(B  •B)u-dS    = 


B^dV 


V 


m 


Taking  the  unperturbed  boundary   (e  s  o)   in  (49)  we  obtain 


(50) 
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li 
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6    P        =      I      -  B    dV    +     (D    ^  (p     -  —  B   )u    dS  (51) 

m  )      M  /     9n    '  o       2/n  n 


V 
m  i 

2 
where  we  have  interpreted  the  operator   u -V   on   p    -B    /2/Li    as  a  normal  deriva- 
tive since    p    -B    /2/u   is  constant  on  the  equilibrium  interface.     The  constraints 
o 

are  not  visible  in  (51)  but  are  implicit  in  the  class  of  vector  fields    B   which  are 

to  be  inserted  in  the  formula. 

2 
The  first  term  in    6"P       is  always  positive,    hence  stabilizing.     The 

?^ 
second  term  is  stabilizing  if  B    /2/li  -  p     increases  going  into  the  vacuum  do- 
main (the  positive  normal  in  (51)  is  directed  outwards}.     This  result  may  also 
be  stated  in  terms  of  0;     we  set 

Q(S)     =    f  (p     -fB^)     .    .A(1b2,   -pM     .  (52) 

9n      o       2/Li  9n    2/Lt  9n 

2 
Since  the  fluid  variation    6    P     is  non-negative,    we  have  stability  if  Q   is  posi- 
tive on  the  whole  interface.     If  Q   is  anywhere  negative,    the  question  seems 
open  to  doubt  because  the  other  terms  compete.     However,    the  vitally  signifi- 
cant fact  about  the  second  variation  is  that  the  sign  of  Q(S)    alone  determines 
the  question  of  stability.     This  follows  from  the  fact  that  one  can  always  find 
a  localized  special  variation,    i.  e.  ,    a  special  deformation  of  the  interface, 

u   (S),    confined  to  an  arbitrarily  small  part  of  the  surface  and  zero  elsewhere, 
n 

which  assigns  a  value  to  the  first  term  of  (51)  smaller  than  the  absolute  magni- 
tude of  the  second  term.     Consequently,    if  Q(St    is  anywhere  negative,    a  varia- 

2 
tion  exists  which  makes    6    P       negative.     Furthermore,    the  special  variation 

can  be  chosen  to  make    1  u   dS  =  0,     so  that  an  incompressible  fluid  variation 

2 
can  be  fitted  to  it  and,    hence,     6    P    ^^  0. 
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Summarizing,    we  have 

THEOREM  1:      A  sufficient  condition  for  stability,   by  the  criterion 
of  the  second  variation,    is  that  Q(S)  >  0   for  all   S.      A  necessary  condition  is 
Q(S)  >  0    for  all   S;    i.  e,   a  sufficient  condition  for  instability  is   Q(S)  <  0   any- 
where on  S. 

The  explicit  construction  of  the  special  variation  mentioned  above  is 

too  lengthy  to  be  given  here.     However,    it  is  easily  described  intuitively.     In 

some  instances  there  exists  a  special  variation  in  which  B   is  unchanged  at 

.2 


Fig. 3 


every  point  and  hence    \B    dV  =  0.      This  can  be  accomplished  by  a  perturba- 
tion of  the  interface  which  follows  a  flux  tube  of  the  original  magnetic  field 

(hence  leaves   B    «  0)   and  leaves  the  volume  un- 
changed (see  Fig.    3).      This  possibility  makes  use 
of  the  fact  that  the  interface  is  known  to  be  analytic 
and  B    can  be  continued  across  the  surface.      In 
such  a  perturbation  all  terms  except  Q(S)   are 
zero.     Such  a  variation  does  not  always  exist.    The 
perturbation  must  follow  a  flux  tube,   and  to  specify 
it  at  one  point  of  a  magnetic  line  determines  it  for 
the  entire  line.     In  an  axially  symmetric  geometry  where  the  lines  extend  in- 
finitely in  both  directions  the  construction  is  possible.     On  a  closed  surface 
such  as  a  torus  the  variation  cannot  be  extended  indefinitely  along  a  magnetic 
line  lonless  the  geometry  is  so  special  that  the  lines  close.     However,    we  may 

take  such  a  variation  and  "taper"  it  gradually  to  zero  along  the  length  of  the 

2 
line.     This  will  not  make  the  first  term  of   6    P       zero,   but  will  make  it  arbi- 

m 

trarily  small  compared  to  the  second  term. 

The  local  nature  of  this  stability  criterion  is  a  property  of  the  free- 
boundary  equilibrium  alone,    (see  Section  10). 
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6.     STABLE  AND  UNSTABLE  CONFIGURATIONS 

The  second  variation  condition  for  stability  has  a  simple  geometrical 

1     2 
interpretation.     For  a  vacuum  magnetic  field,    S/i—B   )  ~  (B-V  B.      Consequent- 

ly,    at  a  point  on  the  interface, 

|-(^B^)     -    n-(B-V:B    -    XB^ 
dn     d. 

where    X  is  the  curvature  of  the  magnetic  line  in  the  interface  which  passes 
through  the  point  in  question.     In  the  case  0=0,     we  note  that  the  principal 
normal  to  a  magnetic  line  on  the  interface  is  normal  to  the  surface,    hence  the 
magnetic  line  is  a  geodesic.     The  sign  of  X.  is  positive  if  the  center  of  curva- 
ture of  the  magnetic  line  is  toward  the  fluid,    and  negative  otherwise.     When 
0=0    we  can  restate  Theorem  1  as 

THEOREM  2:        When   0  -  0,    a  sufficient  condition  for  stability  (posi- 
tive second  variation)  is  that  the  center  of  curvature  of  the  magnetic  lines  on 
the  interface  lies  always  toward  the  vacuum;     sufficient  for  instability  is  a  re- 
verse curvature  anywhere. 

An  immediate  consequence  of  Theorem  2  is 

THEOREM  3:        When   0=0,     there  exists  no  bounded  stable  fluid 
configuration  separated  from  a  magnetic  field  by  an  everywhere  smooth  inter- 
face. 

The  proof  is  simple.     Consider  a  "support  plane"  which  is  tangential 
to  the  boundary  and  leaves  the  entire  fluid  to  one  side.     The  magnetic  line 
which  passes  through  the  point  of  tangency  has  its  center  of  curvature  on  the 
fluid  side  assuring  instability. 
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This  theorem  is  relevant  only  for  a  fluid  domain  bounded  by  one  or 
more  topological  toruses.     A  non-singular  tangential  direction  field  can  exist 

on  no  other  kind  of  surface,    and  a  singularity  of  the  tangential  magnetic  field 

2 
would  violate  the  boundary  condition   B     »  2Aip     »  constant. 

The  term   p90/9n   in  (52)  is  the  normal  component  of  the  external 
force  (per  unit  volume)  at  the  interface.     If  we  allow  an  external  force  field 
to  act  under  the  conditions  of  Theorem  3,   then  stability  can  be  obtained  only 
if  this  force  is  directed  inward  at  every  point  where  the  tangent  plane  is  a 
plane  of  support.     However,    if  one  could  find  such  an  inward  directed  force, 
there  would  be  no  need  for  the  magnetic  field  as  a  containing  mechanism. 

The  curvature  criterion  shows  the  inherent  instability  of  the  pinch 
geometry,    in  which  a  plasma  cylinder  and  surrounding  field  are  enveloped  in 
a  concentric  cylindrical  conductor  (Fig.    4a).     The  curvature  of  the  magnetic 


Fig.  4a 


Fig. 4b 


Fig   4c 


lines  is  clearly  in  the  unstable  direction  for  any  helical  or  azimuthal  field. 
The  only  exception  is  the  limiting  case  of  a  purely  axial  field  (Fig.    4b)  for 
which  the  curvature  is  zero.     The  axial  configuration  is  neutral  since  there 
do  exist  perturbations  not  affecting  the  magnetic  potential.     Such  a  situation 
is  unstable  by  the  criterion  that  large  deviations  from  equilibrium  can  arise 
in  the  course  of  time.     However,    such  neutral  perturbations  must  be  cylin- 
drical,   i.  e. ,    infinite  in  extent  with  straight  magnetic  lines.     For  any  finite 
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perturbation,   the  magnetic  potential  can  be  shown  to  increase.     We  therefore 
consider  this  configuration  to  be  stable. 

The  inversion  of  the  pinch  configuration  consisting  of  a  cylindrical 
tunnel  in  the  fluid  created  by  a  current- carrying  wire  is  stable  (Fig.   4c).     Evi- 
dently the  hole  created  within  a  conducting  fluid  by  a  wire  of  almost  any  shape 
will  be  stable.     In  the  presence  of  a  gravity  potential   0   the  equilibrium  shape 
of  the  hole  will  be  altered.     The  stability  criterion  is  now  that  the  radius  of 
curvature    R  of  the  magnetic  line  at  the  top  of  the  hole  be  less  than  twice  the 
depth   h  below  the  surface.     If  the  wire  is  brought  too  close  to  the  surface,   the 

top  of  the  hole  becomes  un- 
stable and  the  configuration 
presumably  shifts  to  the  one 
of  Fig.    5b. 


If  accelerations  are 
present,    we  would  expect  the 
interpretation  of  the  local 
and  instantaneous  acceleration  as  an  external  force  field  to  give  qualitatively 
correct  insights.     For  example,    in  a  radially  oscillating  pinched  discharge  one 
would  expect  an  inward  acceleration  to  be  destabilizing,    an  outward  one  stabi- 
lizing.    Roughly,   the  effect  would  be  destabilizing  while  the  plasma  is  largest 
and  most  stabilizing  when  the  compression  is  greatest. 

Interesting  configurations  arise  if  more  than  one  wire  is  inserted  in 
the  fluid.     Two  neighboring  parallel  wires  carrying  opposed  currents  are  de- 
picted in  Fig.    6a.     A  very  interesting  configuration  arises  from  four  such 
wires,    alternating  in  current  direction  (Fig.    6b).     A  separated  cusped  config- 
uration appears  in  the  center.     By  removing  the  surrounding  unconnected  fluid 
we  obtain  a  stable  configuration  of  finite  extent  (two-dimensionally).     This  does 
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not  contradict  Theorem  3  because  the  boundary  is  not  smooth;    no  tangent 

plane  is  a  plane  of  support.     From 
the  direction  of  the  magnetic  field 
it  is  clear  that  the  "destabilizing 
curvature"  is  not  merely  concen- 
trated at  the  cusps.     This  and  simi- 
lar two-dimensional  configurations 
can  be  found  by  conformal  mapping; 

a  special  case  is  the  familiar  hypo- 

2/3         2/3 
cycloid  X  +  y  «  1.      A  diverse  family  of  related  two-  and  three-dimen- 

sional stable  cusped  geometries  is  treated  in  an  accompanying  paper. 


ig.6b 


7.     MAGNETIC  VARIATION  IN  THE  LARGE 

From  Theorem  1,    one  would  venture  to  guess  that  the  total  potential 

2 
p  X  p      +  P     is  a  minimum  in  the  large  if  the  quantity  B    /2^  -  p      is  positive 

in  the  entire  vacuum  region;    the  quantity  is  zero  on  the  interface  and  the 

second  variation  condition  is  that  it  increase  toward  the  vacuum.     We  shall 

now  present  a  precise  formulation  of  this  condition.     Here,    for  simplicity, 

we  treat  only  the  case   B     *  0. 

n 

Let   V  »  V      +  V.  denote  the  entire  fixed  domain.     An  admissible  con- 
m  f  

figuration  is  any  domain   V     given  with  the  entire  fluid  state  p(x),    p(x),    etc., 
where  the  fluid  state  is  arbitrary  except  that  it  is  assumed  to  be  compatible 
with  the  original  equilibrium:    the  same  total  mass,    a  total  entropy  which  is 
no  smaller,    or,    for  an  incompressible  fluid,   the  same  volume.     The  magnetic 
field  corresponding  to  this  configuration  is  uniquely  defined  by  the  domain  V 

Berkowitz,  ^J. ,    K.   O.    Friedrichs,   H.    Goertzel,    H,    Grad,   J.   Killeen,   and  E. 
Rubin,    Cusped  Geometries  (in  these  proceedings). 
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and  the  magnetic  constraints   B     =0  and  (9).     We  further  restrict  admissible 
configurations  to  deformations  which  do  not  change  the  topology  of  the  equili- 
brium configuration       Within  this  class  we  further  restrict  our  attention  to  the 
accessible  configurations  defined  below.     Remembering  that    p     is  extended 
into  the  vacuum  domain  by  (25),    we  state 

THEOREM  4       For  an  equilibrium  configuration  with  the  property 

2  ,  — 

that   B    /2/L(  -  p     is  positive  in  some  domain   V   which  encloses  the  equilibrium 

fluid  domain   V.,    the  equilibrium  potential    P     =P      +P^isa  minimum  com- 
I  m        I 

pared  to  the  potential  of  any  accessible  configuration  for  which  the  fluid  do- 
main  V     is  contained  within   V.      If  on  the  contrary,    there  exists  a  domain 
contiguous  to  V      for  which   B    /  2(j.  -  p      is  negative,    then   P      is  not  a  mini- 
mum. 

The  intuitive  significance  of  the  Theorem  is  exemplified  by  the  cusp 

2 
geometry  (Fig.    6b).     The  value  of  B"/2jLi    rises  from  the  value   p      on  the  cen- 
tral fluid  domain,    reaches  a  maximum,    and  drops  back  to   p      at  the  outer 
fluid  boundary.     The  question  of  stability  with  respect  to  finite  displacements 
IS  very  similar  to  that  for  the  container  of  water  depicted  in  Fig,    2b. 

We  shall  only  outline  the  proof.     The  second  part  of  the  Theorem 
(necessity'  is  almost  identical  to  that  for  the  like  part  of  Theorem  2.     The 
first  part  (sufficiency)  is  inherently  more  difficult,    since,    as  is  clear  m  the 
cusped  example     there  may  be  states  neighboring  the  equilibrium  for  which 

a(B  /2ju  -  p  ;/an  <  0. 

o 

o  ' 

We  effect  the  comparison  between  V      and  V     by  introducing  inter- 

o  ' 

mediate  stages.     Let   V      denote  that  part  of  V      not  covered  by   V       and  V 

that  part  of  V     not  covered  by  V    : 

V^    =    V'   nV°    ;  V       =    V°  nv'  (53) 

+  m         f  -  m         f 
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Let   S      denote  that  part  of  the  equilibrium  interface  in  common  with  V    .     On 
S   ,    B     =  V  X  IS  a  tangential  direction  field.     We  assume  (later  verify)  that 
X  can  be  extended  continuously  and  differentiably  into  V      in  such  a  way  that, 
setting  X   =«  V  X  we  have 


X    /2/j     -     p      <    0    ,       in   V, 


+ 


(54) 


The  definition  of  X   is  extended  into  V       by  setting  X  =  B    .      Since 

2  r         o  2     "" 

[X    /2/u-p   ]dV  <  0   and  [(B    ;    /2ax  -  p   ]dV  >  0,    we  have 

v+  °       -  V-  °       ~ 


J 


[^(B°)^  -  p   ]dV    > 
2/j  o  — 


V 


li^'-poi''^ 


tn 


m 


(55) 


If  this  construction  of  X   is  possible,    the  configuration  V„  +  V       is  said  to  be 

f  m 

accessible.     This  inequality  is  valid  for  accessible  configurations.     It  is  con- 
jectured that  m  three  dimensions  all  admissible  configurations  with  V     in   V 
are  accessible.     The  conjecture  has  been  completely  proved  only  for  two- 
dimensional  or  axially  symmetric  equilibrium  configurations.     In  three  dimen- 
sions,  a  certain  domain  between  V     and  V  has  been  shown  to  be  accessible. 

Let  a   denote  an  index  of  range   1,  .  .  .  ,  m;     |3   an  index   m+1,  .  .  .  ,  n. 

The  constant  constraints  (9)  are  the  currents    I     and  fluxes   (Po.      From  (6), 

a  -*-p 

(12),    (29)  we  obtain  (55)  in  the  form 


2iu 


x^^v  -x:«E 


n^^ 


< 


V 


m 


J 

V 


p  dV  -  P 
o  m 


m 


(56) 


The  magnetic  fields   B      and  B'    satisfy  the  same  constraints    I     =  l'       (j)^  "  i.o- 


a 


By  its  construction,    X   has  the  same  constant  current  constraints    I     as    B°, 


a 


fluxes  are  not  defined  for  X   since  it  satisfies    curl  X  =  0   but  not   div  X  =  0, 
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Consider  the  functional 


G[X]    = 


^x^dv-l^?;yx] 


V' 


m 


(57) 


where 


I^[X]     = 


■;l 


X-dx 


(58) 


GFXI   is  to  be  defined  on  the  class  of  all  irrotational  vectors  X  on  V       which 
have  the  given  currents   I  ;     it  is  minimized  on  this  class  by  the  vector  B    . 

o  Of 

Noting  that  1°  »  i'       and  that    M*  =  Y^-  In  ?o  '  5^  t;  I    ^   .    we  conclude  for 
the  particular  X    constructed  above. 


f  X^dV 


-T.Y,^l>.   \  i<B')^-v  -  5;, ,3 


V 


V 


m 


m 


M*[B] 


(59) 


Comparison  with  (56;  yields  the  desired  result, 


M*[B']    + 


p   dV     =    P     [B']    >    P     [B    ] 
o  m  —       m 


V' 


m 


(60) 


8.     VARIATION  IN  TERMS  OF  INDUCTANCE 

There  is  one  class  of  equilibrium  configurations  for  which  stability 
analysis  m  the  large  is  particularly  easy  and  even  reducible  to  finite  dimen- 
sional terms.     A  first  restriction  is  that   B     =0,    0=0.      For  simplicity,    we 

n 


1 


See  footnote,    p.    34. 


54 


assume  also  that  the  fluid  is  incompressible.     The  potential   M'''   is  an  explicit 


reduced  to  a  dependence  on  the  finite  number  of  variables,     L      .     Varying  (7) 


function  of  fluxes,    currents,    and  inductance;    the  dependence  on  the  domain  is 
reduced  to  a  depen 
and  (6),    we  obtain 

I       =    y^  (L      I     +  L      I   )  ,                     I       =   /     ,  (X      I      +  X      i    )       (61) 
—  r         ^ — '        rs  s  rs  s  r         ^ rs-*-s         rs— s 

and 

v — ^  .  .  V — ^ 

M    =     >  ^   (I  L      I     +-I  L      I)     =    /($X      $     +4?^      $^  (62) 

^ — ^      r    rs  s       2    r    rs  s  ^ '     — r    rs— s        2— r    rs-^s 

r,  s  r,  s 

Taking  ^^  I   J      +   ^  I   ^       from  (61)  and  combining  with  (62),    we  obtain 

-2_^^IL      I      =2_vi^^      $       =A^^(I?       -    I   ?    )    =    -^^B^.dS      (63) 
^ — '   2    r    rs  s       ^ — '  2— r   rs— s       ^ — ^  2     r— r  r— r  2u 

r,  s  r,  s  J 

where  the  last  equality  derives  from  (9),    (13),    (14).     All  these  expressions 
must  then  equal  the  first  variation  of   M*    under  the  constraints  (9).     If  a 
given  configuration  has  the  property  that  each  component    L        of  the  indue- 
tance  matrix  is  stationary  for  an  arbitrary  incompressible  variation  of  the 
domain,    it  follows  that    M'''   is  stationary  and  the  configuration  is  m  equili- 
brium.    This  is  a  special  type  of  equilibrium  configuration  which  remains  so 

independent  of  the  values  of  the  fluxes  and  currents;    the  boundary  condition, 

2 
B     =  constant,    is  satisfied  in  every  case.     One  such  vacuum  domain  is  the 

region  between  two  concentric  cylinders.     For  a  purely  axial  field,   an 

2 
azimuthal  field,    or  any  helical  combination  of  the  two,    B      is  constant  on  each 

boundary  component.     Each  of  the  entries  of  the  inductance  matrix  is  station- 
ary for  incompressible  variation  of  either  boundary.     Equilibrium  for  such 
domains  is  a  purely  geometrical  property  and  we  shall  find  that  stability  also 
is  purely  geometrical.     We  call  such  equilibria  universal. 
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The  second  variation  of   M'''    for  a  universal  equilibrium  is  simply 
calculated, 

6    M='=     =    -     >       ^I    L       I        .  (64) 

^ — '  2    r     rs    s 

If  L        is  negative  definite,    then  not  only  is  the  equilibrium  stable,    but  it  is 
stable  no  matter  how  the  magnetic  field  is  inserted  into  the  domain.     For  ex- 
ample,   take  the  domain  between  two  concentric  cylinders  with  the  inner  cylin- 
der a  fixed  conductor,    and  the  outer  cylinder  the  fluid  interface  (Fig.    4c). 
This  configuration  is  stable  for  any  pitch  of  the  helical  magnetic  field. 

For  the  analysis  of  stability  in  the  large,    we  introduce  a  partial 

ordering  in  the  space  of  symmetric  matrices.     If  the  matrix   L'     -  L        is  posi- 

rs        rs 

tive  definite  we  write   L       <  L'     .      In  this  context  the  meaning  of  a  maximum 

rs  rs  ^ 

on  a  given  set  of  matrices  is  clear.     Associated  with  the  class  of  admissible 

domains   V      there  is  a  class  of  inductance  matrices   L'       which  is  represented 
m  rs 

by  a  point  set  in  the  n(n+l)/ 2- dimensional  space  of  matrix  elements.     We  prove 
two  theorems: 

THEOREM  5:     If  a  universal  equilibrium  is  stable  for  one  species  of 
constraints,   it  is  stable  for  all. 

THEOREM  6:     A  necessary  and  sufficient  condition  that  a  universal 
equilibrium  be  stable  is  that  the  inductance  be  a  maximum. 

The  first  theorem  is  an  instance  of  the  results  of  Section  7,   but  a 
direct  proof  in  this  special  case  is  much  simpler.     Suppose,    for  example,   that 
stability  holds  for  constant  flux  constraints;    i.  e.    for  all   ^ 

M'(|)     =   Z^^l/^J^    >     M°(5)     =    X^ll^J^  (65) 

r,  s  r, s 
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We  wish  to  prove  stability  under  other  constraints,    i.  e. 

*         ^ — '      2   3-^^  2   a-^a  *         ^ — '      2  /3-^  a-a 

From  (65),    we  have 

M',  (i' ,1')  =  M'(i')  -  y^  r  ^'  >  M°('3)')  -  y^  r  i'   =  M°(i' ,i°) -Y^  i°^'       (66) 

With  X        fixed,     M   (^ ]   is  convex  in  the  variables  ^    .      Convexity  follows  also 
rs  —  —  r 

for    M   ($    ,^p)  "  ^^  I    $       where  the  equilibrium  values    I  ,    ^o   are  fixed; 
this  implies  a  unique  minimum  at   (j)      «  (j)    ,    i.  e. 


—  a      ^a 


With  (661,    this  completes  the  demonstration.     The  converse  is  easily  demon- 
strated by  constructing  the  appropriate  inequality  for  stability  under  the  mixed 
constraints  condition  and  proving  stability  for  constant  fluxes. 

In  Theorem  6,    a  maximum  inductance  yields  stability  for  constant 

current  constraints  directly  from    M'    =  -  M  =  -  >^  —  I  L      I     and,    hence,    for 

^  2    r    rs  s  ' 

all  other  constraints  by  Theorem  5.     Conversely,    stability  for  any  species  of 

constraints  implies  the  same  for  constant  currents,    and  from  the  form  of  M''\ 

that   L        IS  a  maximum, 
rs 

9.     VIBRATIONAL  STABILITY 

Here  we  consider  only  an  irrotational  incompressible  fluid  flow  and 
stipulate  that  all  the  circulations  (which  are  constants  of  the  motion)  are  zero. 
No  external  potential  is  superposed.     This  problem  is  mathematically  elegant 
because  both  the  fluid  velocity  and  magnetic  field  are  harmonic  vectors,    and 
the  analysis  involves  only  the  interface.     The  fluid  equations  are  given  in 
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terms  of  a  velocity  potential  <f), 

A</)  =0    ,        V  =    V</>    ,        1^   +    |v^   +   p/p     =    0 

2  , 
At  the  boundary,    setting  p  =  B    /2/j,     we  have 

M    +     1     2    ^    _5!    .    0  (67) 

at       2  2{jip 

In  (67)  we  take  the  Lagrangian  derivative   d/dt  »  8/9t  +  v-  7   following  the 
fluid  and  linearize  conventionally  to  obtain 


^    .    ±B     .^^v^j^l     =0  (68) 

at         MP    o     at 


n  an  ^^  2mP  J 


where  B     is  the  magnetic  field  at  the  unperturbed  interface  and  the  condition 
(68)  is  assumed  to  hold  on  the  unperturbed  boundary.     Separating  time  in  the 
form   f(x,  t)  =  e       f(x),     we  then  drop  the  bar  and  obtain  in  (68), 


n  an  ^  2/jp  J 


Jtt>    +    —B    ~    +    v^j-^l     =0  (69) 

^       /up    o    at         ^  ^^  i  o..^  I 


We  choose  not  to  write   aB/at   as   iwB   since  it  is  related  to   v     by  a  time- 
independent  boundary  condition.     The  boundary  condition  (69)  expressed  in 

terms  of  4>   and  a</)/an  =  v     presumably  serves  to  determine  the  harmonic 

n 

function  4)(x). 

For  the  variational  formulation,   we  write,    integrating  by  parts 

K    =  |pv^dV    =  ^p(V<^]^dV    =    |P    J     <^v^dS      . 

Vf  ^f  ^f 

There  exists  a  symmetric  positive  definite  integral  operator  relating  (^   and 

V     =  a^/an   on  the  surface 

n 
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and,   hence. 


0(S)     =         LJS,  S')v   (S')dS' 
I  n 


K    »   -p  V  (S)LJS,  S'}v  (S'ldSdS' 

2  I       n  f  n 


(70) 


From  (50)  and  (51)  the  second  variation  of  M'    is 


U''=     «    6    M^ 


1 

9 

-  B    •  B  u    dS   + 

M      o          n              J 

ai 

„ 

)n  i^  2/Li  J  "n 


dS 


(71) 


where  the  fluid  normal  is  taken  positive  and  the  variation   u     may  be  identified 

with  the  velocity  v   .      In  (71)  we  introduce  the  linear  operator, 

n 


B    -B    X 
o 


L     (S,  S')u   (S'/dS 
m  n 


vy 


Again,    it  can  be  shown  that   L       is  symmetric  positive  definite.     We  consider 

m 

K  and  U'''   in  (70)  and  (71)  as  quadratic  functionals  of  the  surface  function  u   (S) 

n 

* 
and  vary  the  ratio   U    /K   or,    equivalently,    introduce  a  Lagrangian  multiplier 

and  take  the  variation 


6(U''  -  XK]     » 


9  •  1  9B 

-^B     -B   +   iu    -^ 
MO  ju     n    9n 


-    >-P<i> 


6u    dS 
n 


Since     1  6u   dS  *  0   it  follows  that    U''  -  XK   is  stationary  if  the  bracket  is  con- 
stant.    This  constant  can  be  absorbed  in  <f>   as  an  additive  time  function  to  make 
to  make  the  bracket  vanish.     This  is  exactly  the  boundary  condition  (69}  pro- 
vided that  .         . 
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u)       «    -X      . 

In  that  event,     if''  =  XK.      The  eigenvalues  of  U''"    with  respect  to  K   as  unit 
form  are  essentially  the  vibration  frequencies.     If  U''~    is  positive  then  the 
lowest  frequency  is  given  by  the  minimum  of  U''"/K.      If   U''    can  be  made  nega- 
tive at  all,   then  it  can  be  shown  that    U*/K   can  be  made  arbitrarily  large  nega- 
tive and  no  minimum  exists.     This  formal  analysis,    of  course,    needs  to  be 
supplemented  by  an  investigation  of  the  spectrum  of  the  operator. 

Some  analysis  of  this  kind  has  been  made  m  the  compressible  case. 

10.     A  STABILITY  CRITERION  WITH  VOLUME  CURRENTS 

One  of  the  mam  features  of  equilibrium  configurations  m  which  a 
sheet  current  separates  a  field-free  fluid  from  a  vacuum  magnetic  field  is  the 
entirely  local  character  of  the  stability  criterion  (Theorem  ij .     The  violation 
anywhere  of  a  certain  inequality  makes  the  whole  system  unstable. 

The  situation  is  quite  different  for  a  configuration  of  mixed  fluid  and 
magnetic  field  with  distributed  volume  current.     In  this  case,    we  shall  show 
that  a  sufficiently  localized  disturbance  is  always  stable.     This  means  that  in- 
stabilities can  develop  only  because  of  cooperative  behavior  involving  an  appre- 
ciable part  of  the  system. 

f  ds 

The  criterion  is  formulated  in  terms  of  the  quantities   L  «    I  -j — r 

dB  dB    2  "^B-lme  1^1 

and    (>  0  }    where   (-r— )      is  the  largest  eigenvalue  of  the  matrix 

dx   —  dx 

- —  •  -r —  ).      The  product  of  L   with  a  representative  value  of    IbJ    represents  a 
9x.      9x.'  ' 

1  J 

weighted  length  of  the   B-lme   where  the  weighting  factor  may  be  thought  of  as 
See  footnote  2,    p.    30. 
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the  volume  per  unit  length  of  a  flux  tube  following  the  line.     Both,    the  quantity 
dB/dx   and  the  current  density   J   contribute  to  the  rate  of  change  of  B.      We 
prove 

THEOREM  7  For  an  incompressible  fluid,  an  equilibrium  config- 
uration is  stable  against  all  disturbances  confined  to  a  region  R  which  con- 
tains no  closed  B-line,  if  J  does  not  vanish  identically  in  any  subregion  of 
positive  volume^    and  if  along  each  B-line  within   R  the  inequality 

L- Max  Uijl   +  dB/dx  j    <    tt  (72) 

is  satisfied. 

The  criterion  is  sufficient  for  stability  but  of  course  stability  may 
well  occur  when  the  criterion  is  violated.     A  similar  theorem  holds  for  com- 
pressible fluids. 

The  inequality  (72)  can  be  interpreted  as  a  requirement  that  the 
weighted  length  of  each   B-line   in  the  disturbance  region  is  to  be  less  than 
the  distance  in  which  the  fractional  change  in   B   is  appreciable. 

We  establish  Theorem  7  by  deriving  positive  definiteness  of  the  sec- 
ond variation  from  (72);    namely  that 


B^  dV    +    ij. 


Bxj    •    udV    >    0  (73) 


R 


for  arbitrary  continuous  perturbations    u   which  vanish  identically  outside   R 
and  which  differ  from  zero  in  at  least  one  point  of  the  region. 

Using  the  electromagnetic  equations  (2)  (omitting  the  restriction   J=0) 
and  the  relations 
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2m  I       \    (B  "  J)  ■  u  dV 
R 


<        \    B^dV    +    M^  (ux  J)^dV 


R 


R 


r7,^T2  \  f  dB  9b1 


u    ds    > 


2 


i.J=l 


s         -    .2       2 

>7  b     -  a         J 

a  a 


u    ds    ,     (74) 


we  obtain  the  general  inequality  which  yields  the  theorem;     namely, 

I     B^dV    +    M      \     (B  X  J)  .  udV 
R  R 


2 

>    ^ 
-     2 


m!jI 


1 


^.|j|.f    |[L(x,y.z)]2 


:m|j  + 


dB    2 


dx 


■!■■ 


dV 


(75) 


The  stability  criterion  of  Theorem  7  can  be  applied  to  any  sort  of  de- 
vice in  which  all  B-lines  begin  and  end  on  conducting  plates.  If  inequality  (72,j 
is  satisfied,   then  the  device  is  stable 

As  an  application  we  treat  the  stability  of  magnetic  mirror  configura- 
tions with  conducting  end-plates,   that  is     such  axially  symmetric  configura- 
tions as  have  flux  tubes  like  those  depicted  in  Fig.    7.     A  class  of  equilibrium 
configurations  having  this  form  is  given  by  the  equations 


T-.  .,.,/,       V     2,       ,  TTZ   _    ,  iTrr , 

B      =    B   (1  +  Xr   )  -  k  cos —  J   ( ) 

X  o  a      o     a 


B      =  -  k  sin —  J,( j)/i 

r  a      1     a 


(76) 


2XkaB 
p  -  ^B^X[R^2+XR^  -  r^2  +  Xr^)] : 

2       O  ITT 


containing  the  parameters  B   ,    X,    k,    a,    R. 


RJ,( )  -  rJ,( )  cos  — 

la  la  a 
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Three  physical  parameters  may  be  varied:     /3    (the  ratio  of  plasma 
energy  to  magnetic  energy),   the  average  mirror  ratio  (the  ratio  of  maximum 

to  minimum  of    JBj    on  a  B-line 
averaged  over  the  machine)  and  the 
aspect  ratio  (the  ratio  of  the  mid- 
plane  diameter  to  axial  length).     We 
have  prepared  a  table  of  stable  con- 
figurations on  the  IBM  704  at  the 
Institute  of  Mathematical  Sciences 
according  to  the  criterion  (72)  ap- 


Fig.7 


plied  to  the  equilibria  (76).     The  table  shows  that,    for  a  fixed  aspect  ratio,    an 
increase  in  the  average  mirror  ratio  will  decrease  the  maximum  /3    for  which 
the  stability  criterion  is  satisfied.     However,    if  the  mirror  ratio  is  held  fixed, 
then  there  is  an  optimum  aspect  ratio  for  achieving  the  highest   3    compatible 


2.2 
Mirror   Ratio 
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with  the  stability  criterion;    these  values  are  plotted  in  Fig.    8  where  the  maxi- 
mum values  of  /3   are  indicated. 


11.     STABILIZATION  WITH  ALTERNATING  FIELDS 

In  this  section  we  consider  the  problem  of  stabilizing  a  pinched  plas- 
ma by  applying  alternating  external  magnetic  fields.     We  suppose  that  the  ex- 
ternal circuits  are  so  arranged  that  at  each  point  of  the  plasma  boundary  sur- 
face the  magnetic  field  is  of  constant  magnitude  and  rotating  with  constant  angu- 
lar velocity.     The  specific  case  we  treat  is  the  simplest  configuration  which 
displays  significant  features  of  the  problem  and  yet  is  susceptible  to  analysis. 
It  is  the  case  of  a  perfectly  conducting  fluid  supported  against  gravity  at  a  hori- 
zontal plane  interface  by  a  vacuum  magnetic  field  parallel  to  the  interface.     If 
the  field  is  constant,   this  equilibrium  position  is  unstable.        When  the  mag- 
netic field  is  uniform  in  space  and  has  constant  magnitude,    we  shall  see  that 

rotation  about  a  vertical  axis  with  constant 
angular  velocity  has  a  stabilizing  effect. 

We  introduce  rectangular  coordinates 
with  the  y-axis  vertical  and  the  horizontal  x,  z- 
plane  representing  the  plane  of  the  interface. 
In  equilibrium,   the  fluid  is  at  rest  in  the  half- 
space   y  >  0;    the  hydrostatic  pressure  due  to 
the  gravity  potential   gy  is  counterbalanced  at 
the  interface  by  the  magnetic  pressure  from  below.     For  the  sake  of  mathema- 
tical simplicity,    we  assume  the  fluid  is  incompressible. 

We  make  the  assumptions  customary  in  this  kind  of  infinitesimal 


Fig.  9 


1 


See  footnote  3,    p.    29. 
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stability  analysis.     We  imagine  the  perturbed  interface   y  «  f(x,  z,  t)    given  at 
time  t  =  0;    the  problem  is  to  determine  the  subsequent  motion. 

In  the  fluid,    the  velocity   u   is  derivable  from  a  velocity  potential 
0{x,  y,  z,  t)   by  the  relation   u  =  grad  0.      At    y  =  +  oo    we  assume  that  both  u  and 
the  perturbation  in  the  pressure   p   vanish.     The  pressure   p  is  determined 
from  Bernoulli's  law 

I?  *    ^^   ^  £  .   gy    =   c,.>  (76, 

At   y  =  -  00    we  assume  the  magnetic  field  is  equal  to  the  unperturbed 

field  B    (t).     The  direction  of  B    (t)    as  a  function  of  time  will  be  prescribed 

°  °  2 

later.     At  the  interface,    we  have  B     =0   and  B    1 2u  =  p.      The  determination 

of  B    can  be  reduced  to  a  scalar  potential  problem  by  observing  that   B    can  be 

written  as   B  =  B    (t)  +  grad  i//,     where   Atp  =  0. 

2 

The  problem  is  linearized  by  neglecting  the   u      term  in  Bernoulli's 

law  and  the  boundary  conditions  are  imposed  not  at  the  actual  interface  but  at 
the  unperturbed  surface   y  =  0.      Recall  that  the  interface   y  =  f(x,  z,  t)    should 
not  be  prescribed,    for  that  would  fix  the  potentials   0   and  i//   at  once. 

Consider  a  component  of  the  Fourier  expansion  of  y  =  f(x,  z,  t),    say 
y  =  r(t)  sinkx.      For  this  component  we  may  use  the  method  of  separation  of 
variables  to  find  explicit  formulas  for   0   and  ip;    these,    when  put  into  the 
linearized  pressure  equation  at  the  mean  interface   y  =  0,    yield  a  differential 
equation  for   r(t); 


r    + 


(By 

— —    cos    0(ti   -  gk  j  r    =    0      ,  (77) 


) 


where   0(t)    is  the  angle  B    (t)   makes  with  the  x-axis.     This  includes  the  static 

o 

case  where   9(t)    is  constant.     For   d(t)  =  7r/2   we  have  the  well-known  result 
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that  the  growth  constant  of  these  waves  is    i  gk  . 

In  the  dynamic  case,    we  take  constant  angular  velocity  6  -  ut.      Equa- 
tion (77;  then  takes  the  form  of  Mathieu's  equation 

d^r/d0       +    (a  +  bcos2e)r    ^    0      ,  (78) 


where 


2   (    2pM 


2,  2  ,  „2    2 

k  N  ,      B    k 


Ski     ,      b   =  4;    ^     .  (79) 


2p/j  J     '  2      2p/u 

To  investigate  stability  for  a  particular  value  of  uj,     one  must  study 
all  values  of  k,    i.  e.    all  points  of  the  curve  (79)  m  the  (a,b)-plane.     The  values 
of  the  parameters    (a,  b)   which  yield  stability  have  been  charted  m  part. 
Meixner  and  Sch^fke      give  a  stability  chart  for    ;aj,  Jbj  <  50:     for  large  values 

of  a  and  b  the  solutions  of  (78)  have  been  studied  in  very  complete  form  by 

2 
Langer.        By  using  their  results  and  methods  we  find  that  the  application  of 

periodic  fields  does  not  produce  stability.     However,    the  rate  at  which  insta- 
bilities grow  IS  controlled  in  the  sense  that  the  growth  constant  has  an  upper 
bound  independent  of  k.      In  contrast  to  the  static  case,    where  for  large   k  the 
growth  constant  is  proportional  to    V~k    and  may  be  arbitrarily  large,    in  the 
dynamic  case  for  a  small  angular  velocity  w,    the  growth  constant  for  large   k 


is  at  most    gV  p/u  /2  I  B    J.      This  may  yet  be  too  large  for  some  practical 


Meixner,    J.    and  F.    W.    Sch^fke,    Mathieusche  Funktionen  und  SphMroid- 
funktionen.    Springer^    Berlin  (1954),    especially  p.    132. 

2 

Langer,    R.    E.  ,    The  Solutions  of  the  Mathieu  Equation  with  a  Complex  Vari- 
able and  at  Least  One  Parameter  Large,    Trans.    Amer,    Math.   Soc,  ,    36, 
637-695,    (1934). 
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purposes,    but  the  application  of  alternating  fields  does  produce  a  stabilizing 
tendency  in  the  presence  of  acceleration. 
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HYDROMAGNETIC  EQUILIBRIA  AND  FORCE- FREE  FIELDS 

H.    Grad  and  H.   Rubin 


1.     INTRODUCTION.     ELEMENTARY  PROPERTIES 

The  equations  governing  the  equilibrium  of  a  perfectly  conducting 
fluid  in  the  presence  of  a  magnetic  field  are  taken  to  be 

2 
'  grad  p  *  J  X  B  f  grad  (p  +  B    /2m)  «  (B  •  V/B/^i 


I 


curl  B  *  ju  J  or  "S  (11 

div  B  =  0  L  div  B  =  0 


where   p  is  the  fluid  pressure,    B   the  magnetic  field,    and   J  the  current 
density.     These  equations  admit  a  large  variety  of  solutions,    i.  e.    of  equili- 
brium configurations  in  which  a  conducting  fluid  is  balanced  by  a  magnetic 
field.     It  is  our  purpose  to  survey  these  possibilities  with  the  expectation  that, 
if  a  solution  has  been  found  to  exist  mathematically  (and  is,   in  addition,    stable) 
it  can  actually  be  constructed  by  sufficient  exercise  of  experimental  ingenuity. 
In  this  context,   it  is  extremely  important  to  discover  exactly  what  data  should 
be  specified  in  order  to  determine  a  solution  uniquely. 

In  addition  to  certain  general  properties  of  these  equations,    we  shall 
consider  solution  in  terms  of  arbitrary  functions,   the  solution  of  well- posed 
boundary  value  problems,    and  several  alternative  formulations  of  the  problem 
in  terms  of  the  calculus  of  variations.     In  principle,    either  the  differential  equa- 
tions or  the  variational  characterization  can  be  taken  as  the  definition  of  equili- 
brium,   the  two  are  only  approximately  equivalent.     Moreover,    a  certain  type 


The  results  (and  conjectures)  outlined  here  were  presented  in  part  at  various 
Sherwood  meetings  starting  with  one  at  Princeton,    Oct,    1954  (WASH  184,    p.  144). 
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of  variational  formulation  (slightly  different  from  that  treated  in  this  paper) 
can  be  instrumental  in  a  stability  analysis,    which  is  essential  to  give  physical 
meaning  to  an  equilibrium  configuration, 

2 

This  problem  has  been  studied  extensively  in  the  literature,       this 

article  gives  a  survey  of  the  various  possibilities  and  extends  and  general- 
izes most  known  results. 

By  inspection  of  (1),    since  B    and   J  are  perpendicular  to    grad  p.    we 
see  that    p   is  constant  on  B-lines  and  on  J-lines;     equivalently.    the  B-lines 
and  J-lines  lie  on  constant    p   surfaces. 

In  the  case  of  a  unidirectional  magnetic  field  (e,  g.  ,    B    (x,  y))   the  gen- 

2  ^ 

eral  solution  is    p  +  B    1 2fj.  ==  constant.      The  fluid  pressure  is  balanced  by  the 

magnetic  pressure,     either   p   or   B    can  be  given  arbitrarily  and  the  other  one 

"filled  m". 

The  integral  form  jf  (li  is 


(p  +  B"/2/Li)n  -   -  BB     f    dS    =   0  .       d)      B    dS    =^    0  .  (2) 

ju  n  /  -1  „ 


These  equations  are  more  fundamental  than  the  differential  equations  (1).     They 


For  example  see.    Berkowitz,    J.      H.    Grad  and  H.    Rubin,    Problems  m 
Magnetohydrodynamic  Stability  (in  these  proceedings),    Bernstein.    I,    B.  ,    E.    A. 
Frieman,    M.    D.    Kruskal  and  R.    M.    Kulsrud,    An  Energy  Principle  for  Hydro- 
magnetic  Stability  Problems,    Report  No.    NYO-7315,    Matterhorn  Project, 
Princeton  University,    Princeton,    N.    J. 

Some  selected  references  are:     Lust,    R.    and  A..    Schluter,    Z  Astrophys.  ,    34, 
263  (1954);     Truesdell,    C,    The  Kinematics  of  Vorticity,    Indiana  Univ.    Press, 
(1954),     Chandrasekhar,    S.  ,    Proc.    Nat,    Acad.    Sci.  ,    42,    1  (1956);     Chandra- 
sekhar,    S.    and  P     C.    Kendall,    Astrophys.    J.,   J^.    457(1957);     Jorgens,    K. , 
Z  f.    Naturforschg.    13,    493  (1958). 
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are  equivalent  to  (1)  when  the  functions  are  sinooth.     In  addition,    at  a  discon- 
tinuity surface,    the  integral  relations  (2)inipiy  that   B      vanishes  and  that 

2  ^ 

p  +  B    / 2^   IS  continuous.     A  discontinuity  surface  must  be  a  flux  surface; 

otherwise  the  Maxwell  magnetic  stress  tensor  would  not  be  compatible  with 

a  scalar  pressure. 


2.     A  FLUID- DYNAMICAL  ANALOGUE 

Using  an  appropriate  identification  of  variables,    equations  (1)  become 
identical  to  the  equations  of  steady  incompressible  inviscid  flow.     Setting  the 
fluid  density  equal  to  unity,   these  equations  are 

(u  •  V)u    +    grad  p*     =    0 

(3) 
div  u    =    0 


w 


here   p*    denotes  the  fluid  pressure.     The  identifications  are   u  ~B/  •fju   and 


2 
-  p*  ~  p  +  B    /2/u,     the  negative  of  the  pressure   p  in  the  magnetic  case  then 

2 
corresponds  to  the  Bernoulli  constant    p*  +  u    /2.      The  interesting  case  is  ro- 
tational  flow,    for  which  the  vorticity  curl   u  j*^  0    corresponds    to      J  5^  0. 

It  is  interesting  to  consider  a  special  case,    viz.    the  analogue  of  the 
fluid-dynamical  free  boundary  or  cavitation  problem  in  which  an  irrotational 
flow  is  separated  at  a  discontinuity  surface  (vortex  sheet)  from  stagnant  fluid 
or  a  cavity.     The  separation  surface  is  determined  by  the  extra  boundary  con- 
dition   |uj    =  constant.      This  is  mathematically  equivalent  to  a  vacuum  mag- 
netic field   (J  =  0)    separated  at  an  interface  (current  sheet)  from  a  field-free 
conducting  fluid.     Since   p  is  constant  in  the  conducting  fluid,    we  obtain  the 
free  boundary  condition    Jb|    =»  constant. 

Although  the  two  problems  are  mathematically  identical,    it  is  per- 
fectly possible  for  a  significant  fluid- dynamical  problem  to  be  uninteresting 
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in  the  magnetic  analogue  (e.  g. ,    a  jet  from  an  orifice.    Fig,    la)  and  vice  versa 
(e.g.,   the  cusped  equilibrium,      -  Fig.    lb). 


Fig.  la 


»       \   Fig. lb 


3.     SOLUTION  OF  BOUNDARY  VALUE  PROBLEMS 


The  characteristics  of  a  system  of  partial  differential  equations 
give  immediate  qualitative  and  possibly  quantitative  information  concerning 
relevant  initial  value  or  boundary  value  problems.     An  elementary  calcula- 
tion yields  four  characteristics  for  the  system  (1).     Two  of  the  characteris- 
tics are  pure  imaginary  as  for  the  potential  equation.     Corresponding  to  these 
two  characteristics,    one  would  expect  to  prescribe  a  single  scalar  boundary 
value  on  the  entire  boundary  of  the  domain,    e.  g.  ,   the  normal  component   B 
of  the  magnetic  field. 


n 


Berkowitz,   J.,    K.    O.    Friedrichs,    H.    Goertzel,   H.    Grad,    J.    Killeen  and 
E.    Rubin,    Cusped  Geometries  (in  these  proceedings). 

2 

See  Courant,    R.   and  D.    Hilbert,    Methods  of  Mathematical  Physics,    First 

English  edition,    Interscience,    New  York-London,   Volume  II,   to  appear  (in 
German,    Springer,    Berlin,    1937). 
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The  remaining  two  characteinstics  are  real,    viz.    the  B-lines  counted 
twice.     Corresponding  to  each  real  characteristic,    one  would  expect  to  be  able 


S^'  Bn>  0 
^3^  B,  =   0 


Fig. 2 


Fig. 3 


to  specify  a  single  scalar  quantity  at  one  end  of  each  B-line.     In  a  geometry 

as  in  Figure  2  in  which  every  B-line  intersects  each  end  of  the  tube,    one  would 

expect  to  give  both  additional  quantities  at  one  end  or  one  at  each  end.     For 

example,    one  may  conjecture  that  the  following  specification  of  boundary  values 

(in  addition  to  B      on  the  whole  surface)  would  be  appropriate: 
n 


Problem  A 


r 


p  given  on  S 

p  given  on  S     (coinpatibly) 


Problem  A^ 


p  given  on  S 


1 


J     given  on  S, 
n  ^  1 


Problem  A, 


Problem  A 


p  given  on  S 


1 


J     given  on  S^ 
n  ^  ^ 


B         (as  well  as  B   )  given  on 
tan  n 

S     (the  vector  B  is  given) 


In  these  problems,    we  prescribe   p   on  boundary  surfaces  in  such  a  way  that 
the  p-lines  are  "simple"  (not  closed)  as  shown  in  Figure  3.     Moreover,    in 
Problem  A   ,    since   p  is  constant  on  B-lines,    matching  magnetic  flux  on  both 
ends  imposes  a  simple  compatibility  requirement  on   p.     For  example,    if  p  is 
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« 


given  on  S      and  the  lines  of  constant    p  are  specified  on   S        this  condition 
1  ^ 

fixes  the  values  of  p  on  these  given  lines. 

An  iteration  scheme  for  Problems  A     and  A     offering  likelihood  of 
convergence  is  as  follows:    suppose  that,    at  a  certain  stage  of  the  iterations, 
v/e  have  a  vector  field  B   satisfying  div  B  =  0    and  the  boundary  condition  for 
B    .      We  find   p   everywhere  in  the  domain  by  carrying  the  boundary  values  of 
p  along  these  B-lines.     The  component   J      perpendicular  to   B   is  then  obtained 
in  the  domain  from  the  equation   grad  p  =  J  x  B.      We  write   J.j    "  aB   for  the 
parallel  component  and  employ  the  requirement   div  J  =  0   to  obtain  along  each 
B-line  the  ordinary  differential  equation 

B  •  grad  a    +    div  J       =    0  (4) 

for  a.      The  given  "initial"  condition  for   J      at  one  end  allows  this  equation  to 

n 

be  solved  uniquely  on  each  line  so  that    J   is  determined  everywhere.     We  now 

solve  for  a  new  B    from    div  B  =  0,     curl  B  =  J   and  the  boundary  condition  for 

B    ,    and  then  we  continue  the  iterations.  ! 

n  I 

We  next  consider  two-dimensional  problems,    i.  e.  ,    problems  in  which 

j 
no  quantity  depends  on   z.      There  are  several  possibilities      1)  B    can  have  the    i 

single  component   B      and  J  the  two  components    J      and   J   ,     2)    B    can  have  the  i 

z  X  y 

two  components    B      and  B      and   J  the  single  component    J   ;     3)  both  B   and   J 
can  be  general  (three- component)  vectors  depending  on  x   and   y  alone. 

The  general  solution  in  the  first  case  has  already  been  given  explicitly,; 
p  +  B    /2/u  =  constant.      In  the  second  case,    the  number  of  characteristics  is  ; 

three  rather  than  four,    the  B-lines  are  counted  only  once.     The  two-dimen-         , 
sional  analogue  to  problem  A     will  be  considered  later.     Corresponding  to  : 

problems  A       A     and  A     we  have  (see  Figure  4) 

Problem  B  Problem  B    ; 

Ci  ^       O  ■■"■■■"  fr 

p  given  on  C  or  C  Vector  B  given  on  C  or  C 
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In  the  third  (general  two-dimensional)  case^   the  characteristics  are  counted 

as  in  the  full  three-dimensional  case. 
Correspondingly  we  list  (as  above,    the 
analogue  to  problem  A     is  left  to  later): 


Fig. 4 


C,:  Bn<0 
^2-  Bn  >  0 
C3:  Bn  =  0 


Problem  C   : 

p  given  on  F 

J     given  on  F, 
n  ^  1 


Problem  C    : 


p  given  on  F 

J     given  on  F^ 
n  ^  2 


Problem  C 


Vector  B  given  on  F 


The  two-dimensional  problems   C     and  C     are  equivalent. 


With  axial  symmetry,   the  0- coordinate  corresponds  to  the  coordinate 

z   of  the  two-dimensional  case.     We  have  the  same  three  subdivisions:    1)  B   , 

J   ,    J   ;     2)  B    ,    B    ,    J^;    3)  B  and  J  general  (three- component)  vectors  depend- 
r       z  r        z       fl 

ing  on   r   and   z  alon 
alone  which  satisfy 


ing  on   r   and   z  alone.     In  the  first  case,    B„   and  p  must  be  functions  of  r 


(5) 


In  the  second  and  third  cases,    we  identify  problems  D        ,    D     and  E   ,    E   ,    E 
to  correspond  to  B  B   ,    and  C    ,    C    ,    C    . 

Additional  justification  of  the  above  conjectures  will  be  given  later; 
here  they  are  suggested  merely  by  the;  counting  of  characteristics. 

The  problem  of  force-free  fields,    viz.    solution  of  the  system 

curl  Bx  B     «    0  div  B     «    0  (6) 

is  a  special  case  of  the  equilibrium  problem  obtained  by  adopting  the  special 
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boundary  value   p  =  constant.     We  identify  problem  F'    as  the  special  case  of 
A     or  A        F"    as  the  special  case  of  C     or  C    .    and  F'"    as  the  special  case  of 
E     or  E         Problems  F"  and  F'"  and  the  two-dimensional  and  axially  symmet- 
ric  versions  of  F  '. 


4.     CHARACTERIZATION  OF  THE  MAGNETIC  FIELD 
IN  A  PRESSURE  SURFACE 

We  first  introduce  the  concept  of  a  surface  harmonic  vector  field  on 


a  surface   S   (e.  g.  ,    a  constant   p  surface).     A  tangential  vector  field 

vec 
(2) 


(2)  (2) 

X        =    V       0   is  said  to  be  a  surface  gradient    (or  irrotational  vector  field)  if 


f       (2; 
it  is  the  projection  of  a  three-dimensional  gradient  or  if     ©  X        •  dx  =  0    for 

every  closed  curve   C   which  bounds  a  portion  of  S.      The  conjugate  vector 

(2)  (2) 

field  n  X  X        =  Y         (where   n   is  the  normal  to   S)    is  said  to  be  a  surface  curl 

r   —(2) 
(or  solenoidal  vector  field):    we  have      h    Y      ds  =  0   where  v    denotes  the  nor- 

(2\  io\  (2i  (2) 

mal  in   S  to  the  curve   C.      If  X     '   *    V^    '0  =  n  x  V       0   (that  is,    if  X         is  both 

irrotational  and  solenoidal),    we  call  it  harmonic  and  say  that   0   and  i//   are  con- 
jugate surface  harmonics.     In  the  special  case  where  S  is  a  plane,    0   and  \p 
satisfy  the  Cauchy-Riemann  equations. 

(2) 
In  a  simply- connected  plane  domain,    a  harmonic  vector   X     '    is 

uniquely  determined  by  the  boundary  values  of  the  normal  or  tangential  com- 

(2) 
ponents  of  X       .      Specifying  either  at  the  boundary  is  equivalent  (except  for 

a  trivial  additive  constant)  to  specifying  i//   or   0   at  the  boundary.     Exactly  the 

same  facts  hold  on  an  arbitrary  surface   S. 

In  a  multiply- connected  domain,    the  solution  to  such  a  boundary- 
value  problem  is  uniquely  determined  only  when  certain  additional  data  called 

periods  are  prescribed.     These  may  be  specified  values  of  the  circulations' 

(2)  f     (2) 

X       •  dx   on  each  independent  closed  circuit  or  of  the  fluxes       X       ds   on  arcs 

V 
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f 


which  cut  these  circuits  (cf.    Fig.    5).     If  these  periods  are  non-zero,    the  func- 
tions   0   and  ip   are  multiple-valued. 


X^^^dx 


flux  =  Jx^^^ 


ds 


Fig. 5 


Fig,  6 


The  above  theory  can  be  generalized  to  include  weighted  surface 
harmonics  which  satisfy  the  equation 

X  *Y       0»anxY       (//; 


(7) 


here   0   and  i/y   are  conjugate  harmonic  with  respect  to  a  positive  weight  func- 
tion o. 

Since   J   has  no  component  normal  to  a  pressure  surface   S   ,    we  con- 

P 
elude  that   B    is  a  surface  gradient  on  S    ;    i.  e. 

B     «     V^^^0     -    grad0    -    n  |^  (8) 

9n 

In  Appendix  I,    it  i.3  shown  that,    since  div  B  «  0   and   p  is  constant  on  B-lines, 
there  exists  a  function  w   such  that 


B 


(2) 


grad  p  X  grad  d    "    \  grad  pjn  x  V      w 


(9) 


We  see  that  B  is  a  weighted  surface  harmonic  with  weight  J  grad  pj  on  any 
S  .  The  weight  j  grad  pj  arises  in  converting  from  an  actual  three-dimen- 
sional  solenoidal  vector  with   div  B  ■  0   to  a  two-dimensional  area-weighted 
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solenoidal  vector.     We  note  that   B   will  be  an  actual  three-dimensional  gra- 
dient instead  of  a  surface  gradient  only  if  J  =  0   and,    hence,     p  is  constant. 

These  remarks  allow  us  to  extend  the  previously  conjectured  exist- 
ence theorems  to  cases  in  which  the  p-lines  on  the  boundary  are  closed  curves 
(see  Fig.    6).     On  a  surface   S   ,    the  magnetic  field  is  determined  only  when, 
in  addition  to  the  normal  component  of  B   at  each  end,    we  give  a  period,    e.  g. 
the  m.  m.  f. 

T{p)     =    ()    B-dx  (10) 

^  C 

on  a  curve   C    circling  S    .      This  argument  suggests  the  following  modification 
of  problem  A 

Problem  A 

p  given  on  S 
p  given  on  S 
T  given  for  each  p-line 

We  are  now  also  able  to  insert  the  two-dimensional  (C  )  and  axially  symmetric 
(E  )  analogues  of  A  • 

Problem  C  Problem  E 

p  given  on  T  p  given  on  S 

T  given  for  each  p-line  T  given  for  each  p-line 

In  problem  C       t   is  defined  as  the  line  integral  of  B   over  a  finite  distance   z, 
interpreting  the  figure  to  be  periodic  in  the  z- direction.     In  problem  E  , 
T  =  27rB„r.      In  general,    the  value  of  T(p)    can  be  interpreted  as  the  degree  of 
"twist"  of  the  magnetic  lines  on  each  tubular  p-surface.     It  should  be  noted 
that    J      is  an  alternative  way  of  specifying  this  twist. 
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5.     SOLUTION  IN  TERMS  OF  ARBITRARY  FUNCTIONS 

First  we  consider  the  general  two-dimensional  case  in  which  B    ,    B    , 

^  X        y 

B     depend  only  on  x,  y.      Since   div  B  s  o,    we  can  introduce  a  stream  function 

(//{x,  y)    for  the  B   ,     B      components  of  the  field;    we  then  have 

X         y 

B       =    -  ^  B       =    -^  (11) 


X  ay  y       ax 


so  that 


B     =    nx  gradi//    +    nB  (12) 

where  n  =  (0.  0,  1).,      It  can  be  seen  that    p  and  B      are  constant  on  the  curves 

z 

of  constant  ;//   m  the   x,  y   plane.     This  means  that    p  and  B      are  functions 
(possibly  multiple-valued)  of  tp.      The  stream  function  ip{x,y}    satisfies  the  non- 
linear potential  equation 

Ai/y    +    jup'(iA)    +    Mf'(0)     =    0  (13) 

where 

f(0)     =    ^  B^       ;  (14) 

p'    and   f    refer  to  derivatives  with  respect  to  0       In  (13),     p(i//)    and   f(!//)    are 
arbitrary  functions  of  i//       Hence,    our  system  of  four  first-order  equations  has 
been  reduced  to  a  single  second-order  equation  containing  two  arbitrary  func- 
tions.    One  can  expect  a  solution  to  be  determined  by  specifying  the  arbitrary 
functions  and  boundary  values  for  i//   as  for  the  standard  two-dimensional  po- 
tential equation.     For  solutions  depending  on  x   only,    the  equation 

2 

^    +    Aip'((//)     +    Mf'(0)     =    0  •  (15) 

dx 

can  be  integrated  explicitly  giving  many  interesting  equilibrium  configurations. 
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In  the  case  of  axial  symmetry  m  which  B    ,     B        B      depend  only  on 
r,  z,    we  introduce  a  stream  function  (^(r,  z;    such  that 

B       X    i|^  B       =     -^^  (16) 

r  r  9z  z  r  9r 

and,    hence 

B     =    -n'^gradi//    +    nB  (17) 

where  n  «  (n    .n      n   )  =  (0,  1,0).     Here  i//(r,  z!    satisfies  the  equation 

A*i//    +    /nr^p'(0)    +    ^lV{^|J}    «    0  (18) 

where 

2  2 

A*    .   1 i_9_    +    A_  (19) 

^2  r  3r  ^2 

ar  az 

and 

f(i//)     =    -^  r^B^      .  (20) 

2/j  6 

In  the  case  of  cylindrical  symmetry,    in  which  there  is  dependence  on 

r   alone,    we  may  give  any  two  of  p,    i//,    and  B      as  functions  of  r   and  immedi- 

o 

ately  compute  the  other. 

Interesting  special  solutions  of  equations  (15i  and  (18)  can  be  found  by 
separation  of  variables  after  taking  f    and   p'    to  be  linear  in  ip .      In  some 
cases  eigenvalue  problems  arise.     It  is  thus  clear  that  uniqueness  cannot  be 
expected  in  general.     However,    it  is  possible  to  prove  uniqueness  as  well  as 
existence  even  with  quite  arbitrary   f{ip)   and  p(i//~    for  domains  which  are  not 
too  large. 
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setting   p'(ip)  =  0 


The  force-free  special  cases  are  obtained  by  the  simple  expedient  of 
1 


The  reduction  in  terms  of  arbitrary  functions  offers  support  for  the 
conjectured  existence  of  solutions  to  problems  B  C     or  C       D  E     or 

and  F'" 
to  determine  the  functions    piip)   and  f(i//) 


E3.  F" 


In  each  case,  the  boundary  data  over  and  above  B   serve 

•^  n 


6.     VARIATIONAL  ANALYSIS       ADDITIONAL  BOUNDARY 

VALUE  PROBLEMS 


It  IS  well  known  that  solutions  of  the  fluid  free-surface  problem  can 
be  described  as  those  vector  functions   u(xi    which  make  stationary  the  func- 

r  1  2 

tional      —  u   dV   when  the  fluid  volume  is  held  fixed.     Analogously,    in  the  mag- 
J    ^  /^ 

I      2 
netic  case  we  vary     B    /2/Li  dV    holding  the  volume  of  the  magnetic  domain 

fixed.     Or  we  can  drop  the  restriction  to  constant  volume  and^    as  suggested  by 

the  Lagrange  multiplier  rule,    vary 


r 


-^  B^dV    +    p 

2/Ll  o 


r 


r 


dV    -=      1       4-  B^dV    - 
2^ 


p   dV    +    constant 
o 


'V 


■V 


V 


V, 


m 


in 


m 


r 


(21S 


1     2 

( —  B     -  p)dV    +    constant 
2m 


V  +v 

f       in 


Here   V       represents  the  vacuum  (magnetic)  domain  and  V.  the  conducting 
m  I 


In  this  context,    see  also  Furth,    H.  ,    M.   A.    Levine  and  R.    W.    Wanieck,    Pro- 
duction and  Use  of  High  Transient  Magnetic  Fields  II,    28,    11,    Review  of  Sc. 
Instruments  (195?;,    pp.    949-958, 
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fluid  domain,    the  sum   V     +V^  is  fixed.     The  pressure   p  takes  the  value  zero 

m       f 

in  the  vacuum  and  the  constant  value   p      m  the  fluid.     We  shall  see  that,    if  in- 

o 

terpreted  properly,    the  same  functional. 


Q 


^i-' 


p)dV 


(22) 


serves  for  the  general  equilibrium  problem  with  fluid  and  magnetic  field 
mixed,     here   p  as  well  as    B    will  be  a  function  of  x,  y,  z.      The  class  of  ad- 
missible pairs  (B,  p)   allowed  to  compete  in  the  variation  of  the  given  func- 
tional is  restricted  by  the  conditions  that    B   be  solenoidal  and   p  be  constant 
on  the  B-lines, 


I     div  B     =0 

B ■ grad  p    =    0      , 


as  well  as  by  boundary  conditions  that  will  vary  from  problem  to  problem. 


(23) 


1 


A  simple  way  of  incorporating  both  constraints  (23)  is  to  represent 
B   m  the  form 

B     =    grad  p  x  grad  cj  (24) 

(see  Appendix  I).     The  function  w    may  be  multiple- valued  if  the  p-surfaces 
are  not  simple. 

We  perform  the  variation,    obtaining 


WASH  184,    p.    144. 


(|B^-p)dV    = 


V 


(B  -66-  6p)dV 


V 


\  B  •  (grad  6p  x  grad  u  +  grad  p  ^  grad  6u)  -   6pn  dV 


V 

V 


j  6p(J  ■  grad  w-l)  -  6u(J  •  grad  p)  i-  dV 


r 


(25) 


j  6p(B  X  grad  cj)  -   6u}(B  >^  grad  p)  ?  •  n  dS    = 


For  arbitrary  variations  of   6p  and    6u)   in  the  volume  integral,    we  conclude 


J • grad  u    =    1 
J  •  grad  p    =    0 


(26) 


from  which  we  easily  obtain  grad  p  =  Jx  B   as  the  Euler  equation. 

Next  we  turn  to  the  boundary  variation  and  first  consider  the  case  of 

simple  p-surfaces  (Fig.    3),    with   p   given  at  both  ends    S      and  S    .      It  can 

easily  be  verified  that  there  is  no  contribution  to  the  variation  (25)  from   S     on 

which  B    s  0.      On  S      and  S    ,     from    6B     =0   we  conclude  that    6cl)   is  a  function 
n  12  n 

of  p,     however,    we  are  certainly  free  to  fix  the  value  of  u   at  one  end  of  each 
p-line  on   S    ,    and  this  makes    6d  =  0    everywhere. 

We  then  have    THEOREM  1:     In  the  tubular  geometry  of  Fig.    2,    if  p 
is  given  (compatibly)  at  the  ends   S      and  S      m  such  a  way  that  the  p-lines  are 
simple,    then   Q   is  made  stationary  by  any  solution  of  the  system  (1)  which 
satisfies  these  boundary  conditions.     This  theorem  corresponds  to  the  condi- 
tions of  problem  A    . 
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Next  consider  tubular  p-surfaces  (Fig.    6).     Since  w  and  0   (see  equa- 
tion (8))  may  be  multiple- valued,    we  cut  the  domain  on  a  surface   S  which  ex- 
tends from  the  axis  of  the  nested  surfaces    S     to  the  outer  boundary  S     (Fig. 

p  o  _ 

7).     The  boundary  integral  in  (25)  must  now  be  taken  over  both  sides  of  S  as 
well  as   S   ,    S      and  S   .      As  in  the  previous  case,   the  contribution  from   S 


Fig.  7 


Fig. 8 


vanishes.     A  little  manipulation  shows  that  the  contribution  from  the  cut   S  also 

vanishes.     On   S      or  S„,    since   6B     «  0,    we  must  have    6gj  «  f(p).      We  obtain 
1  2  n 

for  this  part  of  the  variation, 
\  f(p)(Bxgrad  p).dS    «       I  f(p)  T(p)dp    «       I     (f2(p)  -  f^(p))  T(p)dp      (27) 


^^h 


S   -S 
2      1 


using  equation  (II.  11)  of  Appendix  II.     The  difference   6u     -  6u     represents  the 

^  J- 

variation  of  the  twist  from   S     to  S     on  a  given  pressure  surface.     In  order  to 
have  Q  stationary  we  must  either  fix  this  twist  beforehand  so  that    6u    -  6w    «0 
or  else  specify  that  T(p)  «  0;    the  latter  is  a  natural  boundary  condition  in  this 
problem.     To  specify  the  twist  for  a  class  of  admissible  vector  fields   B,    we 
require  the  magnetic  lines  which  originate  on  a  ray   C    on  S      to  end  on  another 
given  ray   C"    on   S     (Fig.   8). 

We  have    THEOREM  2:    With  tubular   p  surfaces,     p  given  compatibly 
at  the  ends    S^    and  S^   and  two  given  rays    C    and   C"    identified,    Q   is 
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stationary  if  B   and   p   satisfy  (1).     If  the  twist  is  not  specified,    Q   is  station- 
ary for  a  solution  of  (1)  which  satisfies   T(p)  *  0. 

There  are  really  two  ways  of  specifying  twist,    since  the  specifica- 
tion of  T(p)    or  the  identification  of  two  rays    C    and   C"    are,    m  an  intuitive 
sensCj    equivalent.     It  is  easy  to  alter  the  problem  by  the  use  of  the  Lagrange 
multiplier  rule  so  as  to  be  able  to  specify  T    instead  of  the  twist.     This  varia- 
tional formulation  corresponds  to  problems    A  ,    C     and  E  . 

Some  modifications  of  these  variational  problems  are  necessary  to 
obtain  force-free  fields.     We  drop  the  term   p  m  the  variational  functional, 
and  obtain 

fl       2 
M    *      Ui-  B    dV    .  (28) 

J  2m 

Instead  of  the  representation  B  =  grad  px  grad  w   we  now  take 

B     =    grad  TT  X  grad  u  (29) 

where   -n   is  no  longer  identified  with  the  pressure.     It  is  now  easy  to  verify  the 
following  theorems; 

THEOREM  3:     For  admissible  sets    (B,  tt)    with  simple  7r-surfaces  and 
■n  given  compatibly  on   S     and  S        M  is  stationary  when   B    is  force-free. 

1  ^ 

THEOREM  4:     For  admissible  sets    (B,  tt)    with  tubular    7r-surfaces  and 
TT   given  compatibly  on   S     and  S      as  well  as  the  twist  from   S     to   S        M   is 
stationary  when  B    is  force- free. 

If  one  is  willing  to  place  a  large  burden  on  the  verification  of  the  com- 
patibility of  given  boundary  data,    then  it  is  possible  to  formulate  problems  m 
which   p  [or  tt]   is  given  at  both  ends   S     and  S      without  regard  to  the  simplicity 
of  p- lines.     We  have 
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THEOREM  5:      Consider  admissible  classes    (B,  p)   [or    (B,7r)]   which 
satisfy   div  B  =  0   and   p  [or    tt]    constant  on  B- lines,     p   [or    tt]    given  compati- 
bly at  both  ends,    and  a  fixed  (compatible)  identification  of  the  ends  of  each  B- 
line.     Then  Q  [or    M]    is  rendered  stationary  when   (B,  p)    satisfies     Vp  =  J  >^  B 
[J-^  B  =  0]. 

This  theorem  follows  when  we  note  that  assignment  of  B      and  p  [or 

n 

tt]   at  both  ends  uniquely  fixes  the  correspondence  between  B-lines  if  the  p- 
surfaces[7r- surfaces]  are  simple,    and  if  they  are  not  simple,    the  additional 
specification  of  twist  serves  the  same  purpose.     Physically,    the  boundary  con- 
dition that  the  ends  of  each  B-lme  be  fixed  corresponds  to  a  perfectly  conduct- 
ing fluid  m  contact  with  perfectly  conducting  end  walls    S      and  S    .      In  the 

J.  ^ 

force-free  case,    the  reference  is  to  a  so-called  "pressureless  plasma"  in 
which  the  gas  pressure  is  negligible  compared  to  the  magnetic  pressure,    but 
which  IS,    nevertheless  a  good  conductor,    such  a  fluid  is  physically  realiza- 
ble since  the  conductivity  of  a  plasma  is  independent  of  its  density.     Theorem 
5  can  be  restated  without  mentioning  compatibility  by  giving  B      as  well  as    p 
at  the  end  S      only.     Specifying  the  mapping  of  field  lines  from   S      to   S      then 

-^  J.  ^ 

uniquely  determines   B      and   p  at  the  other  end,     S    . 

It  is  interesting  to  compare  the  variational  problem  for  the  force- 
free  field  with  the  classical  Dirichlet's  principle  which  states  that    M   is  min- 
imized subject  to  given   B      (and   div  B  =  0)    when    curl  B  =  0       The  additional 
requirement  that  the  ends  of  each  B-line  be  fixed  prevents  this  minimum  from 
being  attained  and  yields  curl  B^  B  =  0    as  the  Euler  equation  instead  of 
curl  B  =  0. 

As  a  final  example,    we  take  as  our  domain  the  interior  of  a  topologi- 
cal torus  and  look  for  solutions  which  have  the  outer  boundary  of  the  torus  as 
a  constant  pressure  surface.     No  boundary  values  can  be  given  since  no  B- 
lines  are  accessible.     However,    the  variational  formulation  itself  suggests 
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what  data  are  required  to  obtain  a  well-posed  problem.     As  before,    we  take 
admissible  sets    (B,,  p)    which  satisfy  B  =  grad  p^^gradu),     where  the    p  sur- 
faces now  are  nested  toruses  about  some  closed  curve   C      as  axis.     For  sim- 

o 

plicitv  we  take   p  to  be  monotone,     p    >  p  >  p  ,     where   p      is  the  value  taken 

o  r  o 

on  the  axis    C      and   p     is  the  value  taken  on  the  outer  boundary.     To  make  u 
0*^1  -^     _ 

single  valued,    it  is  necessary  to  introduce  two  surfaces    S      and  S      as  cuts, 
S      is  a  transverse  cut  across  the  torus  (leaving  ends  similar  to   S      and   S 
of  Fig.    2),    and  S„    extends  from  the  axis    C      to  the  outer  boundary.     The  only 
contribution  to  the  variation  (25)  is  on  the  cuts,    we  have 


6Q 


[6p](grad  d  x  B)  •  n  dS 


[6Lj](grad  px  B)  -ndS    .        (30) 


h*h 


h'^h 


Since   p   (therefore    6p)    is  single  valued,     [6p]  -  0.      The  periods   [u]   are  given 

by 


a^(p) 


grad  (J  ■  dx 


'C 


(31) 


where   C     are  the  independent  closed  curves  on  the  torus    S    .      In  Appendix  II 

it  is  shown  that 


o^ip)     =    -    f^   $^(p) 


(32) 


where  ^   (a)    represents  the  magnetic  flux  across  that  part  of  S.    for  which 
p    <  p  <  a.      From  equation  (II.  ID  of  Appendix  II,    we  obtain  the  form 


6Q 


\r^6o^    -    r^6o^\dp 


(33) 
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where   p      and   p,    are  the  values  taken  by    p  on  the  axis    C      and  on  the  outer 
o  1  ■^  o 

boundary  of  the  torus  respectively. 

We  now  state    THEOREM  6:     In  the  class  of  admissible   (B,  p)    defined 
by   div  B  =  0,    B  •  grad  p  -  0    and  given    1   (p),     ?„(p)    (or  a  (p),    a   (p))    for 
p    >  p  >  p  ,    Q   is  stationary  when   grad  p  =  Jx  B. 

We  can  obtain  a  formulation  with  given  t  (p)    rather  than    %_  .(p)    by 
modifying  Q  according  to  the  conventional  Lagrange  multiplier  rule. 

For  force-free  fields,    exactly  the  same  analysis  yields  the  result 
that   Q   is  stationary  if   ^ ,    is  a  given  function  of    ^  9   ^^^^    0       attains  the  fixed 
values    (P  „   and  0    on   C      and  on  the  outer  boundary:    in  parametric  form  we 
give   J   (tt)   and   \  A-n) . 

As  a  special  case  of  a  toroidal  geometry,    we  can  take  a  problem  with 
cylindrical  symmetry  (dependence  on   r   alone)  and  introduce  periodicity  to 
provide  the  toroidal  topology.     In  this  case,    the  above  conjectures  are  trivial- 
ly proved. 

We  summarize  various  sets  of  data  which  are  believed  to  define  def- 
inite equilibrium  configurations.     In  the  tubular  domain  of  Fig.    2,    we  specify 

B      as  shown  and  we  also  specify   p  and   J      at  one  end  or    p  at  one  end  and  J 
n  '^         -^    *^  n  ^  n 

at  the  other.     The  justification  lies  m  the  counting  of  characteristics,    m  a 
heuristically  appealing  iteration  scheme  and  m  a  direct  verification  by  inte- 
grating in  terms  of  arbitrary  functions  in  certain  special  cases.     Or  we  can 
specify   p  at  one  end  and  the  corresponding  terminal  points  of  each  B-lme  m 
a  manner  compatible  with  the  given  values  of  B    .      This  is  confirmed  by  count- 
ing of  characteristics  supplemented  by  the  known  properties  of  surface 


Results  similar  to  these  have  also  been  obtained  by  M.    D.    Kruskal  and  R. 
Kulsrud,    Equilibrium  of  a  Magnetically  Confined  Plasma  in  a  Toroid  (in  these 
proceedings). 
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harmonics,    by  variational  analysis,    and  by  integration  in  terms  of  arbitrary 
functions  m  certain  special  cases.     In  the  toroidal  geometry,    we  conclude  that 
the  two  fluxes    j)    (p)    or  m.  m.  f.  's   T.(p)    can  be  specified  arbitrarily;    the  jus- 
tification is  provided  by  variational  analysis  and  by  explicit  solution  m  special 
cases. 

APPENDIX  I.     REPRESENTATION  OF  AN  INCOMPRESSIBLE  VECTOR  FIELD 

IN  TERMS  OF  TWO  STREAM  FUNCTIONS^ 

We  note,    for  arbitrary  0(x,  y,  z),    i//(x,  y,  z)    and  arbitrary  a{(f),ijj)  (which 
means  that   a   is  constant  on  the  intersection  of  0  =  constant   and  i//  =  constant), 
that 

divlo-  grad  0x  grad  (//)  «  0  (I.  1) 

We  remark  further  that,  if  (1. 1)  holds  for  a(x,  y,  z),  0(x,  y,  z] ,  \p(x,  y,  z)  in  a 
domain  where  0  and  t//  are  independent  functions  (i.  e.  ,  grad  0X  grad  (//  ^  0}, 
then  it  follows  that   a  =  Q'(0,  i//). 

Next  consider  a  solenoidal  vector  field  B   and  a  small  region  of  space 
which  is  simply  covered  by  the  B-lines  which  intersect  a  transverse  surface   S. 
We  have    THEOREM  I.  1      There  exist  functions   0,    <//,    Q'(0.  0)    such  that 
B  »:  a  grad  0  x  grad  i//. 

To  prove  Theorem  I.  1  we  choose  any  independent    0   and  i//   on   S   (i.  e.  , 
0  =  constant   and  <//  =  constant    form  a  coordinate  system;       We  then  carry  the 
values  of  0    and  ip   off  S   on  the  B-lines.     Since   grad  0x  grad  ip   is  parallel  to 
B,     we  conclude  that   B  =  a  grad  0  x  grad  i//   and,    hence,    en  =  a(0^  ip)   by  the  above 
remark. 


The  results  are  well  known,    but  the  details  of  the  proof  illuminate  the  appli- 
cations. 
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We  note  that  the  flux  through   dS   is  given  by 


d(p     =    B    dS    =    d0di//    , 
—  n 


which  justifies  the  terminology  "stream  functions' 


(1.2) 


We  again  take  an  arbitrary  solenoidal  vector  field  B   and  assert 
THEQREMI.  2:There  exist   0,0    such  that   B  =  grad  0  x  grad  i//.      In  fact,    if  0 
is  any  given  function  which  is  constant  on  B-lmes,    then  there  exists  ;//    such 
that   B  =  grad  0x  grad  ijj. 

Choose   0    arbitrarily  on   S   so  that  the  0- curves  are  simple.     Intro- 
duce  s    as  the  arclength  on  a  0-curve,    v    as  the  normal  to  the  0- curve  in   S. 
Note  that  if  B  =  grad  0  X  grad  0   then   B     =  {30  / dv){d4j  / ds) .      This  suggests  that 
we  construct  ip   as 


ip{s) 


1      90 

B      dv 

n 


ds    , 


(1.3) 


integrated  along  each  0-line  (the  value  of  <p   at  one  end  of  each  0-line  can  be 
arbitrarily  assigned).     As  before^    we  carry  the  values  of  0   and  0    off  S   as 
constants  on  B-lines.     From  the  previous  theorem,     B  =a  grad  0Xgradi/y.     By 
the  construction  of  i//,    a  =  1   on   S,     since  a   is  constant  on  each  B-line,    a  =  I 
throughout., 

If  the  0- lines  are  taken  as  closed  curves,    then  the  construction 
yields  a  ip   which  is  not  single- valued,    and  a  cut  should  be  introduced. 

APPENDIX  II.      INTEGRATION  FORMULAS 

We  shall  use  two  basic  formulas:     First,    given  two  functions    f(x'»   and 

g(x)    defined  m  a  space   x  =  (x_  .  .  ,  ,  x   ),    we  have 

I  n 
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f(x'dx    = 


a<g<b 


(11.  1) 


where 


dS 


jrad  g| 


(IL  2) 


The  integral  (II.  1/  over  the  shell   a<g<b   is  written  as  an  iterated  integral 
first  on  the  surface   S   ,    then  with  respect  to   g.        The  second  identity  is, 


grad  0  <  grad  (//    ■dS=      \    dij)  djj    =    (^    0  dip 


i//d0    ; 


(II.  3; 


here  S   is  an  arbitrary  surface  in  three- space  with   C   as  its  boundary.         The 
surface  S    may  have  to  be  cut  to  make   0   and  i//    single  valued.     If  S   is  a  torus, 
on  which  there  are  the  two  independent  closed  curves    C     and   C    ,    this  formula 
reduces  to 


grad  0/  grad(/y    •  dS     =    W^iip]^  -    [(l)]^[ip\ 


(II.  4) 


where 


[0]      . 
1 


d0     , 


[^].     = 


J 


d(// 


(II.  5) 


C 


See  Khinchin,    A.    I.  ,    Mathematical  Foundations  of  Statistical  Mechanics, 
Dover,    New  York,    p.    34. 

2 

See  Blank     A     A,,    K.    O.    Friedrichs  and  H     Grad,    Theory  of  Maxwell's  Equa- 
tions Without  Displacement  Current,    NYO-6486,    Institute  of  Mathematical 
Sciences,    New  York  University  (1957). 
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are  the  periods  of  0   and  i//   respectively. 

We  apply  these  formulas  to  obtain  the  identity 


B   dV    =        \    dp       \   d0dw 
a  <  p  <  b  a  S 


(II.  6) 


for  an  equilibrium   B^     i.  e.    for   B    satisfying 

B    =    V       0     =    grad  p  x  grad  cj 
Using  (II.  2)  and  then  (II.  1)  we  verify 


(II.  7) 


(^^)p 


V       0  •  (grad  px  grad  u)} 


grad  p 


(2) 
grad  u^V       0  ■  dS 


du  d0 


In  the  special  case  of  a  torus,    we  have 


B^dV 


a<p<b 


I  T^(p)CT2(p)    -    r^{p)o^{p)l^  dp 


where 


< 


r 


a.(p)     =    ()      grad  uj  •  dx 
1 


T.(p)     =    a»       B  •  dx 
i 


(II.  8) 


(II.  9) 


(11.10) 
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On  a  transverse  surface   S,    we  now  compute 


i 


(B  X  grad  p- n)dS    «    -      I    T(p)dp 


(11.11) 


a<p<b 


We  have  n  as  the  normal  to   S,    and  introduce  v   as  the  normal  in  S   to  a  p- 
curve,    and  s   as  arclength  along  a  p- curve  (Fig.    9).     The  proof  is 


Fifl.9 


/Bxgradp-n)      »       I    (Bxgradp-n)- 


ds 


o(2) 

V      p 


(11.12) 


(Bx  1/  •  n)ds 


B  •  dx    «    -  T(p)     . 


Again,   on  a  transverse  surface  S, 


(^J 


ds 


»        I  (gradp  X  gradu -n) 
P  1  I  ^<^>p| 


(p      w  •  dx    »    a(p) 
C 


(11.13) 


Consequently, 


$(b)     -     J  (a) 


B    dS    »     I    a(p)dp 


a<p<b 


or 


(11.14) 


CT(p)      »      5  '  (P) 


(11.15) 
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Taking  signs  into  account,   we  can  rewrite  (II.  9)  as 


"a<p<b 
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THERMONUCLEAE.  REACTION  RATE  IN  A  DISCHARGE 

H.    Grad^ 


1,     INTRODUCTION 

The  electrical  conductivity  of  a  plasma  is  found  by  computing  the  de- 
viation from  Maxwellian  of  the  electron  distribution  produced  by  an  externally 

applied  electric  field.     This  problem  has  been  solved  for  the  Boltzmann  equa- 

2 
tion  using  the  Chapman- Enskog  theory  by  Landshoff      and  for  the  Fokker- 

3 

Planck  equation  by  Spitzer  et  al.         In  principle,    the  same  methods  could  be 

used  to  compute  the  ion  distribution  and  from  it  the  thermonuclear  output. 

However,    for  this  application,    the  high-velocity  tail  end  of  the  ion  distribu- 

4 
tion  IS  crucial,    and  this  cannot  be  evaluated  by  the  conventional  methods. 

The  reason  lies  m  the  method  of  linearization  employed.     The  electric  field, 

E,    is  assumed  to  be  small,    and  the  resulting  deviation  of  the  distribution   f 

(0) 
from  the  Maxwellian   f         is  also  small       The  "forcing"  term    E  -  9f/9|    is 


The  author  wishes  to  acknowledge  his  appreciation  to  Albert  A.    Blank  for  his 
valuable  assistance  in  preparing  this  paper. 

2 

Landshoff,    R.  ,    Transport  Phenomena  in  a  Completely  Ionized  Gas  m  Pres- 
ence of  a  Magnetic  Field,    Physical  Review,    1_6,    904,    (1949). 

3 
Cohen,    R.    S.  ,    L.    Spitzer  and  P.    Routly,    The  Electrical  Conductivity  of  an 

Ionized  Gas,    Physical  Review,    80,    230,    (1950),    Spitzer,    L.  .    and  R.    H^rm, 

Transport  Phenomena  in  a  Completely  Ionized  Gas,    Physical  Review,    8^,    977, 

(1953). 

4 
The  same  difficulty  occurs  for  runaway  electrons.     The  methods  used  here 

for  the  tail  of  the  ion  distribution  could  also  be  used  to  obtain  analytic  results 

for  the  tail  of  the  electron  distribution. 
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replaced  by 

d^    _         jE_£    (0) 
91  "     RT 

(see  equation  (12)).     This  is  a  small  perturbation  of  a  Maxwellian  only  if   E  •  ? 
IS  small,    not  when    E  alone  is  small.     Furthermore,    the  collision  term,    which 
is  a  quadratic  functional  of  f  m  either  the  Boltzmann  or  the  Fokker- Planck 
equation,    is  usually  linearized  according  to  the  scheme,     f  =  f        +  f      , 

QU.t)    -Qh'"*,/"')    .    Q(f'»',f"»)    .    Q(f"»,f<">) 

assuming  that   f        «  f  We  shall  find  the  perturbation   f        to  be  much 

larger  than  f  m  the  tail  of  the  distribution  even  though  f  is  m  an  abso- 
lute sense  always  small.  This  fact  vitiates  the  use  of  this  simple  lineariza- 
tion technique. 

We  shall  be  able  to  avoid  these  difficulties  by  making  use  of  certain 
special  features  of  the  Fokker- Planck  equation  and  by  introducing  an  asymp- 
totic expansion  m  the  electron-ion  mass  ratio       In  particular,    the  assumption 
f       «  f         will  not  be  made  m  the  tail  of  the  distribution. 


2.  THE  FOKKER- PLANCK  EQUATION 

The  two  equations    (i  =  +  and  -)    can  be  written  in  the  form 

i  i  2 

1      r    8§  9§    9§        2     rs  8§         r 

m  r  r      s  r 


For  one  derivation,    see  Rosenbluth.    M     N.  .    W,    M     MacDonald  and  D.    L. 
Judd,    Fokker-Planck  Equation  for  an  Inverse-Square  Force,    Physical  Review, 
107,    1,    (1957).     Another  derivation  is  given  m  Grad,    H.  ,    Therinonuclear  Re- 
action Rates  m  an  Electrical  Discharge,    NYO-7977,    Institute  of  Mathematical 
Sciences,    New  York  University,    January,    1958. 
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The  summation  convention  is  used  for  tensor  subscripts,    time  and  space  vari- 
ation have  been  dropped;     f  (?)    is  the  number  density  in  velocity  space; 
e^  =  +  e   where   e   is  the  electronic  charge.     The  dynamical  friction  coefficient 

a      and  dispersion  coefficient   b        are  given  by 
r rs  '^  -^ 


a  (C) 
r 


z 


■      mJ  ^ 


(2) 


b'     (C)     =    l'   >_  B''     (?) 

rs  ^-. — '      rs 

J 


(3) 


where 


T  1  ^  1        ^ 

L       = —    log  A 

47rX^m')^ 


(4) 


X       =       127T 


(KkT) 


3/2 


2n  e 


(5) 


A-^(C) 
r 


V 


V 


I  f'(n^dr7 


(V    =n    -  I 

r         r         r 


(6) 


B^     (?) 

rs 


V"6       -V   V 
rs 


^-^  f^iwdn 


Y" 


(7) 


Rationalized  MKS  units  are  used;     X  is  the  permittivity  of  free  space,    X   is  the 
ratio  of  the  Debye  length  to  the  mean  distance  of  closest  approach  in  a  Cou- 
lomb encounter,    n   is  the  ion  (or  electron)  density.     Other  convenient  forms 
for  the  collision  term  are 


_a_  ,^    i    _af_ 

a?        2     rs  9? 


-  Q^  f^) 

r 


(8) 
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2     rs  8C    a? 
r      s 


r  ac         a? 

r  r 


f^ 


(9) 


where 


r  ^-: — '  i  r 


(10) 


Q       =     L    /        :  A''      . 


(11) 


The  coefficients   A   and  B    can  be  evaluated  explicitly  for  a  Max- 


wellian   f, 


JO) 


n 


(27rRT) 


W2    ^^P^-  2RT^ 


(12) 


In  terms  of 


X 


C      ' 


c   =  Vrt 


(13) 


we  have 


(0)    _       /JjL   _£ 

'r        ~    'Vtt   RT      2 
'  X 


F(x) 


B 


(0) 


rs 


XX  6  3x   X 

/ .  r    s       ,    ,  ,    rs  r   s  ,  ^,    ,  J 

^(6^^--^)e(x>     -     <-r--T-^^^^M    ' 


(14) 


where 


e(xi     = 


x      1     2 
e  dy      , 


0 


X         1     2  12 

_  if       -2^  "2^ 

F(x^     =    -  e  dy    -    e 

0 


(15) 
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A         is  directed  radially  inward.      B         is  not  isotropic;    it  can  be  written 
r  rs 


rs  12  2     rs 


X    X 

r   s 


(16) 


B.    measures  the  dispersion  in  the  radial  direction,   and  B^   measures  the  tan- 
gential dispersion  on  a  sphere  x  *  constant. 

The  dependence  of    |  A       |,    B 
and  B      on  x  is  shown  in  Fig.    1.     We  have 


Fig.  I 


r 


A 


(0) 


[2    n     F(x) 
V  TT  RT     x 


[2       n      2F(x) 


'RT      X 


^       tI'^    -/RT 


(17) 


For  small  x  we  have  the  expansions 


r 


< 


A<0)     . 
r 


1    [2    _n 

■    3   >J  TT     R^ 


rtV^-^^^  ^ 


V. 


„  2      2      n       .^  3      2^ 

B    « (1    -    —X      + 

1        3^,ry^  10 


2      2      n       .,  12^ 

B„  »    —    I (1     -    — -X      +    . 


3    i   TT 


/rt 


10 


.) 


(18) 


and  for  large   x. 


A 


(0) 

V 

X 

n       r 
RT     3 

X 

^1  - 

2n      1 
/RTx^ 

^2    ~ 

Vrt    "" 

X 

(19) 
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the  terms  neglected  are,   in  this  case,   exponentially  small. 

The  coefficients   a      and  b        have  a  simple  intuitive  interpretation. 

r  rs 

A  particle  which  is  known  to  have  the  velocity  f      at  time  zero  "diffuses"  in 

the  interval   dt   into  a  normal  distribution  peaked  about   C    +  a   dt    with  second 

r        r 

moments   b      dt.     The  equilibrium  distribution   f        is  the  result  of  a  competi- 
rs 

tion  between  a  directed  frictional  force,    a   ,    tending  to  bring  each  particle  to 

rest  and  a  random  dispersional  "force"   b        directed  "outward".     In  equili- 

rs 

brium,   the  logarithmic  derivative  of  f  (which  is    -  ?    /RT)   is  given  by  the 
ratio  of  Q  to  b      (b     and  b     are  defined  as  in  (17j). 

Many  interesting  properties  can  be  deduced  by  inspection  of  the  for- 

+ 
mulas  for  a     and  b      .      For  example,   the  ion  friction  coefficient,    a   ,    has  a 
r  rs 

local  minimum  for  a  speed  somewhere  between  the  mean  ion  speed  C     and 
the  mean  electron  speed   C   ;    it  is  dominated  by  ion- ion  collisions  for  low  ion 
velocities  and  by  ion- electron  collisions  for  high.     Although  the  transverse 
component  of  the  dispersion,    b        is  ion  dominated  for  all  velocities,   the  radi- 
al  component,    b   ,    becomes  electron  dominated  above  a  certain  speed  between 
C     and  C   .      One  important  conclusion  we  draw  is  that  the  relaxation  time  for 
that  part  of  the  ion  distribution  which  lies  between   C     and   C     is  on  the  order 
of  the  mean  ion- electron  relaxation  time  (which  is  about    7m    /m     slower  than 
the  relaxation  time  for  the  bulk  of  the  ion  distribution}. 


3.     LINEARIZATION  OF  THE  EQUATIONS 

For  the  mam  part  of  the  distribution,    a  linearized  solution  of  the 
form 

f(?)     =    f^°N?:     +    (E-?^(^(?^/2)  (20) 

is  valid.     Let  us  assume,   tentatively,   that  such  a  solution  has  been  obtained 

for  both  the  ion  and  the  electron  distributions.     Since  A   (?)   and  B      (||   are 

r  rs 
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given  as  integrals  over    f,     they  can  be  computed  from  (20>  and  the  results  are 

valid  for  all  values  of  C        In  other  words,    the  functions    A    (|^    and  B      (?)    are 
r  rs 

perturbed  linearly  in    E   for  all  values  of  ?    even  though  this  is  not  true  of  f(?i 

It  is  assumed  (later  verified)  that  the  tail  of  the  distribution  is  exponentially 

small  even  though  it  inay  differ  markedly  frotn    f  If  we  insert  the  known 

functions   A    {^\    B      (?)    into  the  Fokker- Planck  equation,    the  collision  term 
r  rs 

becomes  a  linear  second- order  differential  operator  on   f.      The  complete 


Fokker- Planck  equation  is  now  linear  (one  need  not  tamper  with  the  term 
E  ■  9f/9C),     and  it  is  uniformly  valid  for  all  values  of  ?        Note  that  the  equa- 
tion is  homogeneous      E   does  not  multiply  an  inhomogeneous  term;     it  is  a 
parameter  in  the  differential  equation      Because  of  this,    we  can  no  longer  ex- 
pect the  solution  to  be  linear  m   E   as  given  m  (20). 

We  shall  find  that  the  deviation  of  the  ion  distribution  from  Maxwelli- 

an  is  extremely  small  (except  for  the  tail).     Consequently,    we  can  take  the 

(0) 


Maxwellian  values   A 

the  electron  coefficients     we  write 


(0(  +  + 

B         for  the  ion  contributions  to   a      and  b    . 


For 


A       .    A^^'     -I-    A^^' 
r  r-  r 


rs  rs  rs 


(21) 


where 


J  V 


r  V    6        -  V  V 
T^^l^  TT         _rs_Vs  ..1     2.    , 

^rs     -"    \       ::3 \^<2^    ^^^ 


V 


(22) 


We  are  interested  only  m  that  part  of  the  ion  distribution  for  which  ?  «  C 
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For  the  electron  coefficients,    (22),    we  use  an  expansion  about  small  ?,    ob- 
taining 


(1)  4 

r  3        1     r 


rs  15        1  rs  r   s  s    r 


(23) 


where 


.00 


(/)(x'  dx 


(24) 


The  value  of  0      can  be  found  without  solving  for  the  entire  electron  distribu- 

2 

tion  function,,    (^(C    / 2) ,    by  integrating  the  ion  Fokker-Planck  equation  after 

multiplying  through  by  C    ■      The  result  is 


[a^(?)    +    ^  E   ]f^(?)dC     =    0 
r  m         r 

+ 


(25) 


Using  the  identity 


A'(?(f'(C!d?     -    0 
r 


(26; 


for  A      and  the  approximation  (23'  for   A    ,     we  find 


2_ 

4tt 


em 

~  2 

JL    m 

+     + 


3K"m 


e   logX 


(27) 


Surprisingly,    the  only  fact  that  we  need  to  know  about  the  perturbed  electron 
distribution  m  order  to  compute  the  ion  distribution  is  the  single  parameter 
<f>        which  represents  the  tnean  frictional  force  between  ions  and  electrons. 


101  - 


This  comes  out  of  the  formal  analysis  of  the  next  section,    but  it  can  also  be 
explained  qualitatively  as  follows. 

The  very  slight  deviation  from  Maxwellian  of  the  ion  distribution  is 
also  easily  explained      To  a  first  approximation,    the  ion  distribution  can  be 
considered  to  be  a  6-function  on  the  scale  of  electron  velocities.     The  electric 
field  accelerates  the  ion  stream  relative  to  the  electrons  until  the  mean  fric- 
tional  force  balances  the  applied  field.     However,    the  individual  ions  all  see 
essentially  the  same  value  of  electron  friction,     consequently  for  each  ion, 
the  electric  field  is  balanced  by  the  electron  friction;    to  this  order,    there  is 
no  deMaxwellizmg  effect.     To  the  next  order  of  approximation,    we  consider 
the  electron  friction  as  a  linear  function  of  velocity  over  the  range  of  ion  ve- 
locities instead  of  a  constant.     To  the  same  order,    the  dispersion  is  still  a 
constant,    since  its  slope  is  zero  at  zero  velocity.     It  is  easily  verified  that  a 
linear  variation  of  friction  combined  with  a  constant  dispersion  is  compatible 
with  a  Maxwellian  ion  distribution.     We  conclude  that  only  the  extreme  tail  end 
of  the  ion  distribution,    where  these  approximations  are  not  valid,    will  exhibit 
any  deviation  from  Maxwellian.     It  is  now  also  clear  why  to  a  certain  order  of 
approximation  the  perturbation  of  the  electron  friction  and  dispersion  coeffi- 
cients can  be  computed  from  the  single  parameter  </)        which  selects  the  point 
on  the  electron-ion  friction  and  dispersion  curve  near  which  the  ion  distribu- 
tion lies. 

4.     ASYMPTOTIC  EVALUATION  OF  THE  DISTRIBUTION  FUNCTION 

We  are  interested  in  ion  speeds  m  the  range   C     «  ^  «  C    .      To  make 
this  statement  precise,    we  introduce  the  expansion  parameter   e, 

m               ^  C 

6  +3                                                              ._. 

=    e  ,           -r    =    e        .                                                      (28) 

m  C 
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For  deuterons,     m,/m     ~  3600      C    /C     ~  60.     and  e~  1/4.      Actually,    there 
are  two  small  parameters,     E  and  e      but  they  will  eventually  be  related.     For 
the  ion  coefficients  we  use  the  high  velocity  approximation  (19)  and  for  the 
electron  coefficients  the  low  velocity  approxitnation  (18)  and  (23).     But  first 
we  transform  to  spherical  coordinates  taking   E  as  the  polar  axis, 

E-C     =    Er  cos  0  ,  ?^     =    r^  (29) 

obtaining 

?  9  h         9  h 

lu     9   f       ..    ^         n    9'^f  2    9    f       .2  ^         n,  e.,af 

-b^  --  -  b3Esme^^  -^  — --  +  [-  .  rp^  .  Ecose(p2  -  -  ]- 

8r  2r     dd 

H^cote.Esine(b3-p2.^]i|i 

^^1  „.  2 ^^2       ^^3 


(301 


[3q^  +  r-^  +  Ecose(4rq2  +  r    "^  +  ^7"  ]f     =    0 


where  we  have  set 


b         =    b^  -^    +    b   <6       -   -^1    +    b,(E  C     +  E  I      -  2E-C-^) 
rs  li;2  2rs  2  3rs  sr  2 


(31) 


The  coefficients   b  ,     p  ,     q     are  functions  of  r   alone.     To  the  order  of  approx- 
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imation  taken, 

,  +         nL^    (2  2    /T     3  8    '''^+'^1  ,^  , 

+  X 


+  p2 

+          nL  (   1           2    /T      3  ^     8  ""    +"^1   _  .  )  ,  +           8        ,       + 

b       =    -p-  ?  -  +   ^     -    e  +    -—  (E  x)K  b  -■    --  TT.^    L       , 

2           C  (  x          3  ^  TT                  15        n  )           3          15        1 


(32)  -  continued  on  next  page 
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-nL 


+ 


C 


+ 


1  12      3) 

X 


qt    =    0 


^o 


3^*1^    '' 


(32) 
(continued* 

-  e/m    , 


+     _    J.    /2      3   nL 

^      "     3  J  TT    ^  _3 


C 


+ 


3-^^L    /e 


=    el 


m 


where 


X    =    l/C 


(33) 


Factoring  out  the  Maxwellian, 


JO; 
f    =    f      g      . 


(34) 


we  obtain 

ax  '     X   e      96 

m  terms  of  the  dimensionless  electric  field 


r  3.^   ,3/2   *1 

G    =     -  -(27ri 


n  5  3 

ne   log  A 


(36) 


This  equation  is  to  be  solved  subject  to  the  restrictions    1  «  x  «l/e      and 
£  «  1,      If  X  ~  l/e   or  smaller,    it  is  easy  to  see  that    g   does  not  depart  appre- 
ciably from  unity.     In  this  range,    the  linearized  solution,    (20),    remains  valid, 

2 
In  the  more  interesting  range,    for   x  »  l/e    (e.g.    x  ~  l/e    i,    the  dominant 

terms  m  the  equation  are  those  involving   dg/dx  and   g;    the  solution  is 


2     3    3  r 
exp  ;  -  —  e   X    ^  cos 


X  ~  1  /e' 


(37) 
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The  ditnensionless  electric  field  5    is  essentially  the  energy  gain  of  a  parti- 
cle in  such  a  field  over  a  mean- free- path  compared  to  the  mean  thermal 

energy,     kT.      If  this  quantity  is  large,    no  linear  theory  is  adequate.     The 

3  2 

exponent  m   g  has  order    8,1^   ,    say    1/e"    if  4  ~  e,     and  this  can  be  large, 

producing  an  enormous  perturbation  over  the  Maxvvellian.      On  the  other  hand, 

,4  ,2 

the  exponent  in  the  Maxwellian  itself  has  order    1 /e      for   x  ~  1/e    ,     so  the 

total  distribution  function   f  is  very  small  in  this  range 

One  of  the  most  interesting  features  of  the  formula  (37'>  is  the  minus 
sign  m  the  exponent.     This  implies  that  the  ion  distribution  function  is  de- 
pressed below  the  Maxwellian  in  the  direction  of   E  and  is  enhanced  m  the 
opposite  direction.     Of  course,    since  we  have  already  observed  that  to  low- 
est order  the  forcing  term    E   is  exactly  cancelled  by  the  increase  in  electron 
friction     it  is  no  longer  a  matter  of  simple  intuition  which  direction  this  higher- 
order  effect  will  take.     As  it  turns  out,    it  is  the  fact  that  the  electron  disper- 
sion coefficient  decreases  in  the  direction  of  E  that  produces  the  effect  (the 
ion  distribution  has  been  shifted  m  the  same  direction  as    E  from  the  maxi- 
mum of  b) 

5.     THERMONUCLEAR  OUTPUT 

Since  we  have  computed  only  the  dominant  term  (i,  e.   the  exponent) 
of  the  distribution  function,    we  carry  out  the  thermonuclear  computation  to 
the  satne  order       If  we  write 

1  2 

o    ~     exp;Li(-  mV    »      ,  m  -  m  (38  > 

for  the    D-D   cross-section   (V   is  the  relative  velocity  as  m  equation  (6V),    the 
thermonuclear  output  to  dominant  order  is  given  by 

exp    |M(|mV^:j    f^°  (|*g(§'f^°\€^)g(C^)d?dC^      ..  (39) 
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The  integrand  is     with  the  exception  of  /.( ,     the  exponential  of  a  polynomial  in 

?    and  C         We  introduce  the  variables   V  ■-  c,-?,     W  -■  ?  +C    and  find  that  the 
1  I  1 

(li  11' 

integrand  has  a  maximum  at   W        and  V        given  by 


W 


(1^ 


U(V^°^^^C^ 


ju'(-m(V      )    t 


(40) 


1     r.        1   .2,^  (0).2    6,2 
—  [1  -   -  2=    (V       )    e    /C_^] 


2kT 


.(0) 


V'"     IS  the  maximizing  value  for    i^^  =  0.      We  approximate  the  integral  by  the 
maximum  value  of  the  integrand  to  obtain 


(0)  (.1      ,,,(0>  2^3  ,     2  2  ,( 

V  --  V        exp   '{-      (V       )     e    /CI     > 


(41) 


y 


(0) 


IS  the  Maxwell-averaged  thermonuclear  output,    with    ^  ~  0       It  is  con- 


venient to  convert   6  ■'■  kT   into   kev      n   into   ions /cm      and    E  into   volts /cm 
For  temperatures  which  are  not  too  high,    we  have 


m(V^°*)^/kT    ~     12    4/0'^^ 


(42) 


from  which  we  obtain 


(0)  (  rO    9  <  10^^  Ee^^^,2) 

V        exp   1  L \ ;i J    \ 

(  n  log  X  ) 


(43) 


(44) 


A  representative  set  of  values  for  which  the  exponent  m  (43)  has  order 
unity  might  be  n  ~  10  ions /cm  E  ~  1  volt /cm,  0  ~0  5  kev,  log  X  is  about 
14.  The  most  interesting  point  about  the  formula  for  v  is  the  extremely  sensi- 
tive dependence  on   n,     E      and   6 ,    the  transition  from  a  negligible  correction  to 


1 


Thompson,    W     B    _    Thermonuclear  Reaction  Rates.    Proceedings  of  Physical 
Society,    London,    B70,    1      (1957). 
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an  enormous  one  requires  a  small  change  in  the  parameters.     A  specific  fea- 

2 
ture  is  the  emphasis  on  low  density,     combined  with  the  factor   n      which  is  m 

V       ,    we  find  that   v    has  a  minimum  as  a  function  of  n   when  the  exponent  is 
unity,    and  then  v    rises  very  rapidly.     For  example,    m  a  discharge  of  vary- 
ing density,    it  is  likely  that    E  would  be  approximately  uniform,    the  thermo- 
nuclear output  might  then  be  concentrated  m  the  lower  density  regions.     How- 
ever,   one  should  recall  the  fact  that  the  appropriate  relaxation  time^    i.  e.   the 
time  required  to  build  up  the  tail  of  the  distribution  is.    for  one  thing,    long 
(long  compared  to  the  ion-electron  relaxation  time";,    and  for  another,    inverse- 
ly proportional  to  the  density. 

Formula  (44)  is  also  of  interest  since   W         measures  the  mean  shift 
m  the  center  of  mass  of  the  reacting  deuterons.     The  energy  spectrum  of  the 
neutrons  emitted  has  a  shift  of  the  order  of  the  thermal  energy  when  the  dimen- 
sionless  term  has  order  unity;     m  the  example  above,    the  shift  is  about  0.  5  kev. 
We  recall  that  this  shift  should  be  observed  in  the  direction  opposite  to   E. 
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CUSPED  GEOMETRIES^ 

J.    Berkowitz     K.    O.    Friedrichs     H     Goertzel, 
H.    Grad,    J.    Killeen     E.    Rubin 

1.     INTRODUCTION 

We  are  interested  in  the  containment  and  thermonuclear  possibilities 
of  a  large  family  of  stable  magnetohydrodynamic  free  boundary  equilibrium 
configurations.     The  free  boundary  is  a  mathematical  idealization  in  which 

there  is  a  perfectly  conducting  plasma  containing  no  magnetic  field,    separated 

2 
at  an  interface  (surface  current)  from  a  vacuum  magnetic  field.        The  plasma 

2 
pressure,     p,    is  balanced  by  the  magnetic  field  pressure,    H    /Stt.      Although 

the  original  motivation  for  studying  this  model  was  its  mathematical  simplici- 
ty,  the  sharp  separation  was  soon  found  to  offer  practical  advantages  in  per- 
formance,   provided  that  it  can  be  realized  experimentally.     The  mathematical 
simplicity  allows  some  parts  of  the  theory  of  cusped  geometries  to  be  ana- 
lyzed m  much  greater  detail  than  for  any  other  containment  geometry. 


It  can  be  shown  that  any  finite  plasma  configuration  of  the  free  bound- 

3 
ary  type  bounded  by  a  smooth  interface  is  unstable.        However  there  do  exist 

a  large  number  of  absolutely  stable  configurations  bounded  by  cusped  surfaces. 


The  theory  described  in  this  paper  was  originally  presented  at  Sherwood  meet- 
ings in  Berkeley,    February  1955  (WASH-289,    p.    115),    m  Princeton,    October 
1955  (TID-7503,    pp.    238,    319).     m  Gatlinburg,    June  1956  (TID-7520,    pp.    148, 
373,    394,    400);    in  Berkeley,    February  1957  (TID-7536,    p.    200). 

2 
Grad,    H.   and  H.    Rubin,    Hydromagnetic  Equilibria  and  Force- Free  Fields, 

(in  these  proceedings). 

3 

Berkowitz,    J.  ,    H.    Grad  and  H.    Rubin,    Problems  in  Magnetohydrodynamic 

Stability,    (in  these  proceedings). 

-   108  - 


These  have  been  shown  to  be  stable  for  arbitrary  perturbations,    even  of  finite 
amplitude,    regardless  of  the  equations  which  are  assumed  to  govern  the  plas- 
ma motion.        The  basic  two-dimensional  prototype 
is  shown  in  Fig.    1.     The  magnetic  field  is  a  conse- 
quence of  the  four  line  currents,   alternating  in  di- 
rection,   coupled  with  the  plasma  sheet  currents 
which  take  alternate  directions  on  adjacent  seg- 
ments.    There  is  a  one-parameter  family  of  plas- 
mas of  increasing  size  (Fig.    2)  for  a  given  coil  con- 
figuration.    The  shapes  can  be  computed  by  con- 
formal  mapping.     The  limiting  shape  for  small  dimensions  is  a  hypocycloid, 

2/3         2/3 
X         +  y  *  constant,    and  the  largest  configuration  which  does  not  spill  out 


(a) 


is  given  by  the  equation 


(  (-g-)  cos  ITS  ) 


where  tan  6   is  the  slope  and   s   is  the  arclength.     This  equation  holds  for  the 


1 


See  footnote  3,    p,    108. 
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portion  of  the  curve  in  the  upper  part  of  the  first  quadrant;    the  complete  shape 
is  found  by  symmetry.     The  shapes  can  also  be  computed  nuinerically  for  cur- 
rent arrangements  more  realistic  than  line  currents. 

Among  the  theoretical  questions  that  should  be  asked  about  a  given 
configuration  presumed  to  be  operating  as  a  thermonuclear  device  are  its  sta- 
bility,   rate  of  loss  of  particles  and  of  energy,    and  experimental  feasibility. 
To  answer  all  but  the  first  question,   the  original  idealization  must  be  modified 
to  take  into  account  a  continuous  transition  from  plasma  to  vacuum.     The  pic- 
ture we  have  in  mind  is  still  that  of  a  thin  transition  zone  of  minimum  thickness 
separating  a  reasonable  approximation  to  a  vacuum  from  an  almost  field- free 
plasma,   but  whenever  possible  we  shall  try  to  relate  this  simplified  model  to 
the  more  general  case  of  completely  intermixed  plasma  and  field.     Of  special 
interest  for  comparison  purposes  is  the  opposite  extreme  to  the  free  boundary 
case,    viz.    a  plasma  of  such  low  density  that  the  plasma  currents  have  no  effect 
at  all  on  the  magnetic  field  of  the  external  coils.     One  such  case  is  the  picket 


(b) 


Fig.  3 


fence.       Fig.    3a,    which  has  the  free  boundary  analogue  shown  in  Fig.    3b. 
(This  interface  is  also  explicitly  computable  using  conformal  mapping.) 


'Tuck,    J.    L.,    WASH-184,    p.    77. 
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The  major  advantage  of  the  cusped  configuration  is  stability.        The 
chief  disadvantage,    at  least  of  the  basic  geometry  of  Fig.    1,    is  the  large  rate 
of  loss  of  particles.     One  can  generalize  the  concept  of  magnetic  tnirror  to 
apply  it  to  the  motion  of  a  particle  approaching  a  cusp  even  though  the  magni- 
tude of  the  magnetic  field  does  not  increase  towards  the  cusp,   but  is  a  con- 
stant on  the  interface.     Nevertheless,   because  of  the  extreinely  large  effective 
mirror  ratio  (see  Appendix),    it  is  found  that  the  particle  losses  increase  rap- 
idly with  temperature  instead  of  decreasing  as  do  the  conventional  rnirror 
losses.     For  preliminary  experiments  with  relatively  low  temperature  plas- 
mas,  the  cusps  are  advantageous.     For  eventual  thermonuclear  temperatures, 
a  crossed  field  inserted  in  the  plasma  is  very  effective  in  reducing  losses. 
The  prototype,    shown  in  Fig.    4,    has  a  uniform  longitudinal  field  existing  only 
inside  the  plasma,    and  the  cusp  field  outside. 


Crossed 
field 


Fig.  4 


Fig.  5 


Another  possible  disadvantage  of  a  cusped  geometry  is  uneconomical 
use  of  magnetic  field  if  the  plasma  is  to  be  compressed.     A  possibility  for 
ameliorating  this  situation  is  to  compress  within,    cf.    Fig.    5.      The  tunnel 


This  was  the  first  thermonuclear  containment  geometry  which  was  shown  to 
be  stable,  WASH-289,  p.  115,  and  it  has  been  shown  to  be  stable  under  much 
wider  conditions  than  any  other  geometry  (see  footnote  3,    p.    108). 
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created  by  such  an  internal  conductor  is  stable  against  all  deformations. 
More  conventionally,    field  shaping  by  use  of  appropriately  designed  conductors 
and  by  use  of  iron  in  the  magnetic  circuit  can  be  effective.     This  question  of 
magnetic  efficiency  is  not  entirely  one-sided  since  most  other  configurations 
require  significant  internal  magnetic  fields  in  the  plasma  to  obtain  stability. 
Some  of  these  considerations  have  also  been  studied  by  M.    Levine  (private 
communication). 

2.     SUMMARY  OF  GEOMETRIES 

There  exists  a  large  variety  of  three-dimensional  stable  cusped  con- 
figurations suggested  by  the  prototypes  depicted  in  Figs.    1,    3b,    4,    and  5.     The 
spindle  shape  of  Fig.    6a  is  obtained  by  rotating  the  two-dimensional  shape  of 


Fig. 6 

Fig.    1  around  an  axis  through  two  cusps.     The  resulting  figure  has  a  circular 


See  footnote  3,    p.    108. 
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line  cusp  and  two  point  cusps.     The  coil  arrangement  consists  of  two  opposite- 
ly oriented  Helmholtz  coils.     This  three-dimensional  free  boundary  has  been 
computed  numerically  and  does  not  deviate  much  from  the  rotated  two-dimen- 
sional shape.     By  placing  a  number  of  spindles.    Fig.    6a,    end  to  end,    joining 
and  widening  the  point  cusps,    we  obtain  Fig.    6b.     The  infinite  periodic  limiting 
case  can  also  be  described  as  a  rotation  of  Fig,    3b  about  its  longitudinal  axis. 
Fig,    6c  is  obtained  from  Fig.    6b  by  bending  it  into  a  torus.     The  toroidal  shape. 
Fig.    6d,    is  similarly  obtained  by  rotating  Fig.    1  about  a  horizontal  axis  in  its 
plane;    the  coil  system  consists  of  four  ring  currents.     Fig.    6e  shows  how  a 

crossed  field  might  be  introduced  into 
Figs.    6a  or  6b;    the  field  due  to  the  axial 
wire  is  allowed  to  penetrate  into  the  plas- 
ma.    The  addition  of  an  axial  surface  cur- 
rent (as  in  a  pinch  discharge)  to  Fig,    6b 
modifies  the  shape  sotnewhat;    the  cross- 
section  midway  between  two  circular  cusps  (i.  e.   under  the  external  coil)  be- 
comes narrower.     This  can  become  unstable  for  a  sufficiently  large  axial  cur- 
rent.    Specifically,    since  the  interface  has  negative  curvature,    there  are  two 
asymptotic  directions  separating  directions  of  positive  and  negative  curvature. 

Instability  is  reached  when  the  axial  cur- 
rent rises  sufficiently  to  make  the  total 
current  cross  an  asymptotic  direction, 
see  Fig.    7.     An  axial,    pinch-type  current 
can  be  introduced  in  Fig.    1  (i.  e.   a  cur- 
rent on  the  interface  into  the  plane  of  the 
figure).     The  asymmetric  Fig.    8a  results. 
As  the  axial  current  is  increased,    it  is 
likely  that  a  transition  will  occur  to  the  configuration  shown  in  Fig.    8b. 


Fig.  6e 


stable 


Fig.  7 
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3.     PARTICLE  LOSSES 

Particles  can  be  lost  from  a  system  by  diffusion  across  the  magnetic 
field  and,    in  open-ended  systems,    lost  along  the  magnetic  field  (end  losses). 

The  latter  is  our  concern 
here.     It  is  illuminating 
to  compare  two  extreme 
cases,    collision-  dominated 
mirror  losses  and  non- 
adiabatic- dominated  cusp 
losses.     We  shall  find  that 
there  is  a  continuous  tran- 
sition from  one  to  the  other. 

First  consider  the 
elementary  containment 
argument  for  the  axially  symmetric  geometries.  Figs.  9a  and  9b.  One  as- 
sumes that  the  motion  is 
adiabatic,  i.  e.  the  mag- 
netic moment  /j  is  a  con- 
stant. This  is  the  ratio 
of  the  gyration  energy,  w^, 
belonging  to  the  component 
of  velocity  perpendicular  to 
B,    to  the  magnitude  of  B, 


Fig. 8b 


Fig.  9a 


Fig. 9b 

and  it  is  constant  to  good  approximation  if  the  magnetic  field  seen  by  the  parti- 
cle changes  only  by  a  small  amount  in  one  gyro  oscillation.     The  gyration  ener- 
gy increases  as  the  particle  moves  toward  larger   B,    but  it  cannot  exceed  the 
constant  value,    w,    of  the  total  kinetic  energy.     Consequently,   reflection  occurs 


Post.    R.    F.,   Summary  of  UCRL  Pyrotron  (Mirror  Machine)  Program,    UCRL- 
5044,    (June  1958). 
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when  B    reaches  the  value   w/^.      Furthermore,   the  guiding  center  of  the  parti- 
cle is  constrained  to  lie  on  a  fixed  flux  tube.     This  allows  us  to  conclude  that 
successive  reflections  will  occur  at  exactly  the  same  value  of  z   (Fig.    9).     All 
particles  with   Wj_/w   less  than  a  certain  critical  value  (this  defines  the  loss 
cone)  are  instantly  lost,    and  the  others  are  contained  indefinitely.     Under  the 
assumption  that  the  configuration  is  exactly  axially  symmetric,    there  are  two 
mechanisms  which  can  shift  a  particle  into  the  loss  cone,    collisions  and  alter- 
tion  of  the  magnetic  moment.     The  two  mechanisms  are  independent  and  either 
can  be  dominant;     it  is  even  possible  for  each  one  to  be  dominant  for  different 
classes  of  particles  m  a  single  system.     For  example,    m  Fig.    9b,    an  orbit 
can  be  quite  accurately  adiabatic  for  a  particle  whose  gyro  radius  is  small 
compared  to  its  distance  from  the  axis,    the  loss  mechanism  is  collisions.     On 
the  other  hand,    the  motion  will  be  strongly  nonadiabatic  if  the  orbit  passes  close 
to  the  origin,    and,    for  these  orbits,    collisions  can  be  neglected.     This  effect  can 
also  be  present  m  Fig.    9a  (it  is  always  present  m  Fig.    9b)  if  the  plasma  pres- 
sure is  comparable  to  the  magnetic  pressure  near  the  center  of  the  machine, 
thereby  reducing  the  magnetic  field  m  the  center  to  a  relatively  low  value.    More 
precisely,    we  suppose  that,    on  the  axis  in  Fig.    9a,    the  mirror  ratio  (ratio  of  ' 

the  maximum  to  minimum  value  of  B)   is  large.     A  particle  whose  turning  point       i 
is  close  to  the  position  of  maximum   B    will  have  a  very  small  value  of  Wj_  /w 
when  it  reaches  the  center  of  the  machine.      It  can  move  a  large  axial  distance 
in  one  gyro  period  and  therefore  be  strongly  non-adiabatic.  I 

An  ex..e.e  case  of  nonaa.abat.eUy  .  .He  cu.ped  geo^e.r,.     Howeve.,  J 

we  shall  first  describe  it  in  entirely  different  terms.     We  start  with  the  sim-  1 

plest  model     viz.    with  a  sharp  discontinuity  separating  field-free  plasma  from  I 

vacuum.     At  thermonuclear  temperatures,    the  mean-free- path  is  essentially  '■ 
infinite  compared  to  the  dimensions  of  the  apparatus.     A  particle  trajectory 
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consists  of  straight  lines  joined  by  cycloidal  arcs  (Fig,    10a);    the  limiting  case 
of  zero  gyro  radius  is  the  billiard  ball  model.    Fig.    10b.     In  the  latter  case 


Fig.  10 


(b) 


there  are  no  losses;     every  particle  is  turned  back  before  it  reaches  the  cusp 
unless  it  is  on  the  axis  and  aimed  directly  at  the  cusp.     In  the  former  case,    we 
(tentatively)  define  a  particle  to  be  lost  if  it  reaches  the  cusp.     It  is  found  that 
those  particles  which  are  aimed  close  enough  to  the  cusp  axis  are  lost;    the 
others  are  reflected  back.     But,   it  is  easy  to  see  that  the  aspect  with  which  a 
reflected  particle  approaches  the  next  cusp,   at  the  opposite  end  of  the  machine, 
has  very  little  correlation  with  its  behavior  in  the  first  cusp.     To  good  approx- 
imation,   we  may  assume  that  the  particle  distribution  is  isotropic  in  velocity, 
i.  e.  the  loss  cone  is  refilled  as  fast  as  particles  are  removed  (this  is  in  fact 
a  general  property  of  any  configuration  with  a  very  large  mirror  ratio).     It  is 
possible  to  treat  the  motion  of  a  particle  when  it  is  near  a  cusp  by  studying  a 
suitable  adiabatic  invariant  ju*    instead  of  the  magnetic  moment  /j    (see  Appen- 
dix).    The  "mirror  ratio"  appropriate  to  the  adiabatic  invariant  n*   has  the 
order  of  the  square  root  of  the  ratio  of  the  gyro  radius  to  the  apparatus  dimen- 
sion.    Although  n      is  approximately  constant  while  the  particle  is  near  the 
cusp,   it  loses  almost  all  memory  of  its  former  value  in  crossing  from  one 
visit  with  a  cusp  to  the  next. 
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The  loss  rate  computed  on  the  basis  of  an  isotropic  velocity  distribu- 


1    . 


tion      is  an  upper  estimate.     It  is  a  good  estimate  in  a  strongly  nonadiabatic 
device.     In  other  cases,    the  loss  rate  is  governed  by  the  rate  at  which  parti- 
cles diffuse  into  the  loss  cone  by  collisions  and  by  nonadiabaticity. 

The  precise  computation  of  particle  losses  is  a  very  complex  prob- 
lem involving  the  solution  of  a  self- consistent  field  problem  for  the  transition 
zone  between  the  plasma  and  field,    the  merging  of  two  such  zones  at  the  cusp, 
and  the  interaction  between  these  boundary  layers  and  the  space  charge  sur- 
rounding the  plasma.     We  shall  present  a  sequence  of  models  of  increasing 
complexity  to  estimate  the  loss  rate. 

The  boundary  layer  is  formed  of  two  classes  of  particles,    "bound" 
particles  which  belong  to  the  plasma,    i.  e.    which  follow  orbits  as  shown  m 
Fig.    10a,    and  "free"  particles  which  spiral  nearby  without  entering  the  plas- 
ma proper       The  combination  produces  a  certain  electromagnetic  field  which 
we  must  find,    at  least  to  some  approximation.     It  can  be  shown  that  the  loss 
rate  per  unit  length  of  the  cusp  is  given  by 

uf  du  dv  dw  (1) 


V 


rv   f  du  dv  dw 


u>0 


u>0 


where  f(u,  v,  w)    is  the  particle  velocity  distribution    (u    is  directed  toward  the 
cusp)  and   r(u,  v,  w)    is  the  range  of  a  particle  (See  Fig.    10a),    i.  e.   the  distance 
travelled  in  the  direction  of  the  cusp  in  one  encounter  with  the  boundary  layer. 
The  length   r    can  be  interpreted  as  the  width  of  a  "hole"  near  the  cusp  through 
which  particles  stream  freely. 


Suggested  for  the  picket  fence  geometry.    Fig.    3a,    by  H.    Snyder,    private 
communication. 
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To  start,    we  ignore  the  electric  field  and  assume  that  the  magnetic 

field  jumps  discontinuously.     The  trajectories  can  be  evaluated  explicitly,    and, 

2  2  2  2 

for  an  isotropic  monoenergetic  (i.  e.    fixed  speed   u     +  v     +  w     -  V    ]    distribu- 
tion,   we  compute 
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where  n   is  the  number  density.     At  a  given  temperature,    electrons  and  ions 
are  lost  at  the  same  rate;    the  electrons  move  faster,   but  the  "hole"    r    is 
smaller  by  the  same  factor.     This  picture  must  be  made  self  consistent  with 
respect  to  current  and  with  respect  to  charge.     The  second  is  far  more  im- 
portant.    Its  effects  can  be  approximated  by  recomputing  the  boundary  layer 
orbits  m  the  presence  of  a  constant  electric  field   E   whose  value  is  then  ad- 
justed to  make  the  mean  penetration  of  ions  and  electrons  equal  (the  Debye 
length  is  assumed  to  be  small  compared  to  the  boundary  layer  thickness). 
The  result  is,    approximately. 
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2       2 
For  deuterons,     M/m  =  V    /V     ~  3600,    the  ion  loss  rate  is  cut  down  by  a  fac- 
tor 35  and  the  electron  rate  increases  by  only  a  few  percent.     The  boundary 
layer  is  now  approximately  self-consistent,    but  the  difference  in  loss  rates 
must  be  equalized  by  external  space  charge  effects.     The  precise  mechanism 
is  vital  since  a  possible  factor  35  is  at  stake.     If  the  space  charge  field  does 
not  penetrate  into  the  plasma,    it  merely  accelerates  the  lost  ions  but  returns 
electrons  to  the  plasma,    equalizing  the  net  loss  rate  at  the  lower  value,    v    . 
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The  possibility  of  space  charge  penetration  into  the  plasma  can  only  be  dis- 
cussed realistically  if  the  magnetic  field  is  allowed  to  permeate  the  plasma. 
A  rough  analysis  of  this  problem  yields  the  result  (mainly  in  consequence  of 
the  very  large  mirror  ratio)  that  the  equalization  occurs  almost  entirely  by 
reducing  the  electron  loss  rate  without  appreciably  altering  the  ion  loss  rate. 

A  more  serious  consideration  is  the  possibility  that  some  of  th*e 
electrons  which  are  turned  back  to  the  cusp  may  miss  the  cusp  and  become 
"free"  electrons  as  previously  described.     These  electrons  will  be  trapped 
in  the  external  space  charge  field  and  will  oscillate  between  cusps  in  the  outer 
reaches  of  the  boundary  layer  allowing  the  ion  penetration  to  increase,   there- 
by increasing  r      and  v    .      This  matter  can  be  analyzed  to  some  extent,    and 
arguments  given  to  indicate  that  the  boundary  layer  widens  only  slightly  as  a 
result  of  particle  collisions  and  the  effect  is.  not  large.     However,    we  merely 
summarize  by  stating  that  the  correct  loss  rate  is  somewhere  between  the  two 
values,    v^  and  v   ,    quoted  in  (3),    with  a  fair  chance  that  it  will  turn  out  to  be 
closer  to  the  smaller  value.     An  alternative  description  is  that  the  loss  rate 
is  proportional  to  the  thickness  of  the  boundary  layer  which  thickness  is  some- 
where between  the  electron  and  ion  gyro  radii,   but  may  very  well  be  closer  to 
the  former.     The  boundary  lawyer  thickness  can  possibly  be  controlled  experi- 
mentally,   e.  g.    as  shown  in 
Fig,   11.     An  arrangement  is 
shown  whereby  electrons 
might  be  stripped  from  the 
outer  edges  of  the  boundary 
layer. 

The  Maxwell-aver- 
aged loss  rate,   using  the  more 
optimistic  value  v   ,   of  (3)  is 


g.  M 
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V     = (4) 

640     eB 

From  this  we  find  the  containment  time  of  the  spindle  shape  m  Fig.    6a  to  be 

T    =   10      -Y-  (5) 

where   t    is  in  seconds,     H   is  in  gauss,     R  is  the  radius  in  cm.  ,    and  T   is  the 
temperature  m  electron  volts.     (The  density  is  determined  by  the  temperature 
and  magnetic  pressure  which  equals  the  gas  pressure.  )     For  example,    if 
H  =  10,  000  gauss,     R  =  100  cm.  ,     T  =  10  ev.  .     we  have  a  containment  time  of 
one  second.     In  the  thermonuclear  range,     H  =  100,  000  gauss,     R  =  100  cm.  , 
T  =  10,  000  ev.    yields   t  -  10  milliseconds.     In  present  day  low  temperature 
(<  1  kev.  )  plasmas,    the  containment  time  of  a  cusped  configuration  can  be 
much  longer  than  the  time  one  can  maintain  the  pulsed  magnetic  fields. 

The  general  case  with  field  and  plasma  intermingled  such  as  m  Fig. 
9b  is  extremely  complicated.     The  outer  layers  with  relatively  small  mirror 
ratio  are  collision  dominated  with  a  relatively  low  loss  rate  and  with  electron 
losses  preferred.     There  is  possibly  an  intermediate  layer  m  which  electrons 
are  adiabatic  and  ions  are  not  since  the  latter  have  larger  gyro  radii.     The 
loss  rate  in  this  layer  could  show  a  preference  for  either  ions  or  electrons 
depending  on  the  parameters.     An  inner  layer  might  be  nonadiabatically  dom- 
inated for  both  ions  and  electrons  with  the  electron  losses  predominating 
again. 

The  space  charge  problem  in  the  general  case  is  formidable.     Al- 
though rough  estimates  of  loss  rates  can  still  be  made,    it  would  be  advisable 
to  look  to  experiment  for  guidance. 

Next  we  briefly  consider  the  effect  of  a  crossed  field.    Fig.    4,    and 
first  take  the  case  of  a  sharp  discontinuity  m  the  magnetic  field.     Computations 
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involving  generalized  adiabatic  invariants  such  as  those  described  in  the  Ap- 
pendix have  been  made  to  compute  the  loss  rate  through  a  cusp.        The  im- 
portant point  to  realize  is  that  this  computation  (which  yields  a  lower  loss 
rate  than  without  the  crossed  field)  applies  only  to  the  plasma  layer  which  is 
within  one  gyro  radius  of  the  vacuum.     The  inner  particles  suffer  no  losses 
at  all.     More  precisely,    the  losses  from  the  outer  layer  can  only  be  replaced 
by  diffusion  across  the  magnetic  field.     This  considerably  lower  loss  rate 
then  becomes  the  dominant  one.     In  the  case  of  a  gradual  transition  instead  of 
a  discontinuity  of  the  magnetic  field  the  situation  is,    of  course,    very  complex, 
but  the  general  qualitative  features  can  easily  be  ascertained.     The  mirror 
reflection  argument  can  be  carried  out  largely  as  m  the  case  without  the 
crossed  field.     Moreover,    the  magnetic  field  never  vanishes,    and  nonadia- 
batic  effects  are  considerably  reduced. 

We  can  summarize  the  effect  of  the  crossed  field  as  follows.     In  the 
general  case  it  reduces  nonadiabatic  losses  to  the  generally  lower  level  of 
coUisional  losses,    and  m  the  sharp  boundary  case  to  the  still  lower  level  of 
diffusional  losses  across  the  magnetic  field. 

4.     THERMONUCLEAR  POSSIBILITIES 

For  simplicity  we  consider  only  the  geometry  of  Fig.    6a.     The  state 
of  the  plasma  is  given  by  three  parameters,    e.  g.    volume,    pressure,    and 
temperature.     Actually  we  use  the  radius.     R,     of  the  circular  cusp  section  m 
cm.  ,    the  magnetic  field,     H,     in  gauss     and  the  temperature      T,     m  electron 
volts.     The  energy  loss  through  the  cusp  m  watts  is  then  (from  equation  (4)}, 

P        =    0    08  RTH  (6) 


Killeen,    J.,    TID-7536,    p.    200. 
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One  can  estimate  the  Ohmic  losses  as,, 

ri 

and  the  thermonuclear  output  (energy  of  charged  particles  only)  for  a  deuteri- 
um-tritium mixture  at  20,  000  ev.    as 

P^     =    lO'^^R^H^     .  (8) 

T 

Clearly,     P^    can  be  made  to  dominate   P„  +  P_,   by  making   R   and   H   large 
(the  ratio  of  bremstrahllung  losses  to  thermonuclear  power  depends  on  tem- 
perature alone). 

Let  us  define  "breakeven"  by  the  condition    P      =  3(P     +P    ).      A  pos- 

5  9 

sible  set  of  breakeven  figures  is    H  =  10    ,     R  =  40,     P      =  6  x  10    .      The  density 

is  about  10       ions/cm    .     The  power  output  would  be  comparable  in  size  to  a 
conventional  large  generating  station.     We  envision  a  cycle  somewhat  as  fol- 
lows.    Let    <S    represent  the  original  energy  of  the  plasma.     At  some  later 
time  (m  this  example,    about  one  millisecond)  when  6    has  been  reduced  by 
one-half  due  to  cusp  losses,    the  thermonuclear  input  has  reached   3£'/2   bring- 
ing the  total  plasma  energy  to   2^  and  the  plasma  pressure  to  double  its  orig- 
inal value.     An  adiabatic  expansion  against  the  magnetic  field  then  transfers 
the  amount    I-   into  the  generator  which  supplies  the  field.     The  remaining 
energy  (including  the  considerably  greater  energy  carried  away  by  the  neu- 
trons) could  be  removed  by  a  heat  cycle  of  lower  efficiency.     Of  course,    these 
numerical  results  should  be  considered  only  as  representing  the  correct  order 
of  magnitude. 


Post,    R     F.  ,    Controlled  Fusion  Research  -  An  Application  of  the  Physics  of 
High  Temperature  Plasmas,    Reviews  of  Modern  Physics,  ^^    pp.    338-362 
(1956). 
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A  number  of  points  should  be  mentioned  here.     One  property  of  cusp 
losses  is  that,    in  a  given  device,    high  speed  particles  are  lost  more  quickly, 
tending  to  deplete  the  Maxwell  taiL     In  particular,    this  applies  to  the  charged 
reaction  products.     On  the  other  hand,   the  hot  reaction  products.,    m  modify- 
ing the  boundary  layer,,    serve  to  reduce  the  loss  rate  of  the  original  "cold" 
plasma  and,    in  particular,,    of  its  Maxwell  tail..     The  use  of  the  optimistic  loss 
rate  formula  (4'  is  probably  much  more  than  offset  by  the  gain  that  would  arise 
from  the  use  of  crossed  fields. 

Essentially  all  of  the  methods  of  initiating  and  heating  a  plasma  can 
be  used  m  a  cusped  geometry.      Four  parallel  discharges  alternating  m  direc- 
tion might  be  set  up  in  a  longitudinal  magnetic  field  to  initiate  Fig.    4;    or  the 
coils  in  any  of  the  cusped  configurations  could  be  pulsed  after  preionization  or 
after  creation  of  the  plasma  by  a  Shockwave:     or  the  plasma  could  be  inserted 
in  an  already  existing  vacuum  magnetic  field  of  appropriate  type  by  an  arrange- 
ment of  plasma  guns.     A  particularly  intriguing  possibility  is  to  pulse  the 
magnetic  field  in  one  shot  with  enough  energy  to  bring  the  plasma  to  thermo- 
nuclear temperature.     The  plasma  will  undergo  dainped  oscillations       Although 
the  initial  acceleration  would  probably  be  too  high  to  maintain  stability,       the 
acceleration  would  very  likely  drop  to  a  stable  value  within  a  cycle  or  so. 
Other  confinement  geometries     e.  g.    the  stabilized  B    -pinch,    are  not  suffi- 
ciently stable  to  allow  such  violent  treatment       Roughly  speaking,    in  exchange 
for  a  high  degree  of  stability  we  compromise  somewhat  m  particle  losses, 

APPENDIX       THE  CUSP  AS  A  MIRROR 

The  theory  of  adiabatic  invariants  of  a  Hamiltonian  which  is  slowly 
varying  in  time  can  be  modified  so  as  to  apply  to  a  Hamiltonian  which  has  an 


See  footnote  3,    p.    108. 
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"approximately  ignorable"  coordinate.         This  allows  the  containment  argu- 
ments of  section  3  to  be  extended  to  an  axially  symmetric  system  of  the  type 
of  Fig.    9a  even  when  the  gyro  radius  is  comparable  to  the  radial  dimensions 
provided  only  that  it  is  small  compared  to  the  axial  length   (z   is  the  approxi- 
mately  ignorable  coordinate).     In  this  case  the  adiabatic  invariant,    m    »    is 
w, /cj     where  u*   is  the  frequency  of  oscillation  of  the  particle  in  the  "frozen" 
or  cylindrically  symmetric  Hamiltonian  (with  a  truly  ignorable  coordinate 
corresponding  to  the  instantaneous  value  of  z)   instead  of  the  magnetic  mo- 
ment which  is   w, /u.     (The  gyro  frequency,    w,    equals   B    except  for  the  con- 
stant factor   e/m.) 

Exactly  the  same  remarks  apply  to  two-dimensional  problems.  For 
a  particle  approaching  a  line  cusp,  the  frozen  Hamiltonian  corresponds  to  the 
problem  of  two  parallel  plates.    Fig.    12.     With  no  electric  field,    u*   has  the 
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Fig.  12 


d  is  the  separation  between 

the  plates  (which  varies  slow- 

1        2 
ly  with   z)   and  vv     =  —  mv.    . 

With  the  adiabatic  invariant 
/i*  =  w.  /u*,    we  can  analyze 
the  motion  of  a  particle  in  either  a  cusp  or  a  conventional  mirror.     The  mir- 
ror ratio  would  be  the  ratio  of  u*   at  the  cusp  (viz.    u)   to  its  value  elsewhere; 
for  large  d,    this  ratio  is  essentially  the  ratio  of  the  gyro  radius  to   d/n.     A 
mirror  ratio  can  be  assigned  to  the  cusp  as  a  whole  by  estimating  at  what 
value  d  the  particle  makes  its  last  encounter  with  the  magnetic  field  before 
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crossing  over  to  the  next  cusp.     This  mirror  ratio  turns  out  to  be  on  the  order      1 

of  M  =  1/2  -fa.    where   a   is  the  ratio  of  the  mean  range,    r,    to  the  length  of  the       ! 

i 
figure.     The  corresponding  correction  to  the  loss  computation,    i.  e.   the  meas- 
ure of  anisotropy  of  the  velocity  distribution,    is  the  small  number --j-j:(1  +  log 8 M). 
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